clear all
close all

dt=0.2*1e-6;

Tstop=50e-6;

C=1e-6;
L=1e-6;
w0=1/sqrt(L*C);
R0=sqrt(L/C);
vcm=10;
ilm=vcm;

A=[0 ,-1/C; 1/L,0];

Tfe=eye(2)+dt*A;
Tbe=inv(eye(2)-dt*A);
Ttr=inv(eye(2)-(dt/2)*A)*(eye(2)+(dt/2)*A);

nstep=ceil(Tstop/dt);

t=dt*(0:nstep);
phi=w0*t;
vc=vcm*cos(phi);
il=ilm*sin(phi);
tplot=t/1e-6;

xfe=zeros(2,nstep+1);
xbe=zeros(2,nstep+1);
xtr=zeros(2,nstep+1);

xfe(1,1)=10; % vc(0)
xfe(2,1)=0; % il(0)
xbe(1,1)=10; % vc(0)
xbe(2,1)=0; % il(0)
xtr(1,1)=10; % vc(0)
xtr(2,1)=0; % il(0)

for i=1:nstep
    xfe(:,i+1)=Tfe*xfe(:,i);
    xbe(:,i+1)=Tbe*xbe(:,i);
    xtr(:,i+1)=Ttr*xtr(:,i);
end

vcfe=xfe(1,:);
ilfe=xfe(2,:);
vcbe=xbe(1,:);
ilbe=xbe(2,:);
vctr=xtr(1,:);
iltr=xtr(2,:);

figure(1)
title('forward Euler')
subplot(2,1,1)
plot(tplot,vcfe,'b',tplot,vc,'r')
axis([0, 50, -30,30])
xlabel('t [us]')
ylabel('v_C [V]')
subplot(2,1,2)
plot(tplot,ilfe,'b',tplot,il,'r')
axis([0, 50, -30,30])
xlabel('t [us]')
ylabel('i_L [V]')

figure(2)
title('backward Euler')
subplot(2,1,1)
plot(tplot,vcbe,'b',tplot,vc,'r')
axis([0, 50, -30,30])
xlabel('t [us]')
ylabel('v_C [V]')
subplot(2,1,2)
plot(tplot,ilbe,'b',tplot,il,'r')
axis([0, 50, -30,30])
xlabel('t [us]')
ylabel('i_L [V]')

figure(3)
title('trapezoidal')
subplot(2,1,1)
plot(tplot,vctr,'b',tplot,vc,'r')
axis([0, 50, -30,30])
xlabel('t [us]')
ylabel('v_C [V]')
subplot(2,1,2)
plot(tplot,iltr,'b',tplot,il,'r')
axis([0, 50, -30,30])
xlabel('t [us]')
ylabel('i_L [V]')

%          1      2     3     4     5    6      7    8   9
outdata=[tplot' vcfe' ilfe' vcbe' ilbe' vctr' iltr' vc' il'];
save outdata outdata