clc;//очистили экран командой clc. Далее вводим исходные данные
E=200; R1=10; R2=40; R3=40;  R4=10;L=0.4;C=30/1000000; uCo=200;iLo=0;
p=poly(0,"p"); //Объявили переменную р аргументом полинома
//Далее рассчитали коэффициенты матрицы контурных сопротивлений
Z11=R1+R2+1/(C*p);Z22=R1+R3+R4;Z33=R3+L*p;
Z12=R1;Z21=R1;Z13=0;Z31=0;Z23=-R3;Z32=-R3;
//Далее объявили матрицы контурных сопротивлений и контурных ЭДС
Z00=[Z11 Z12 Z13;Z21 Z22 Z23; Z31 Z32 Z33];E11=uCo/p;E22=E/p;E33=0;
Ek=[E11;E22;E33]; Ikt=Z00\Ek // Рассчитали матрицу контурных токов
IcP=Ikt(1,1),ILP=Ikt(3,1),//Присвоили  значения 1 и 3 элементов
//Далее объявили iC1 -  числитель, iC2 - знаменатель
iC1=numer(IcP);iC2=denom(IcP);
P00=roots(iC2),Nn=size(P00,'r'),//Нашли корни знаменателя.Nn-число корней
iC3=derivat(iC2); // Нашли производную от знаменателя
iC1o=coeff(iC1),iC2o=coeff(iC3);
//Далее методом Хевисайда находим начальное значение 1 экспоненты
k=[1:1:Nn];
iCt_1=sum(iC1o(k).*(P00(1))^(k-1))/sum(iC2o(k).*(P00(1))^(k-1))
iCt_2=sum(iC1o(k).*(P00(2))^(k-1))/sum(iC2o(k).*(P00(2))^(k-1))
t=[0:0.001:0.1];
ioC=iCt_1*%e^(P00(1)*t)+iCt_2*%e^(P00(2)*t);
plot(t,ioC,'b');xgrid();