function x=heun(A,B,fw,fw2,dt,t,x0)
maxit=t/dt;
x=x0;
[m,n]=size(A);
I=eye(m,m);
f1=inline(fw);
f2=inline(fw2);
for i=1:maxit
   u1=feval(f1,i*dt)
   u2=feval(f2,i*dt)
   u=[u1;u2]
   p=(I+dt*A)*x+dt*B*u;
   x=x+0.5*dt*(A*x+B*u+A*p+B*u);
   iter(i)=i;
   zmx1(i)=x(1);
   zmx2(i)=x(2);
   zmx3(i)=x(3);
end;
plot(iter,zmx1,iter,zmx2,iter,zmx3)

% R1=100
% R2=100
%C1=100e-6
%C2=220e-6
%L=1e-3
%x0=[0;0;0]
%E=[R1*C1 0 L; 0 R2*C2 L; C1 C2 0]
%A1=[-1 0 0;0 -1 0; 0 0 1]
%B1=[1  0;0 1;0 0]
%A=inv(E)*A1
%B=inv(E)*B1
%fw='10*sin(2*pi*50*t)'
%fw2='20*sin(2*pi*50*t)'
%dt=1e-5
%t=1e-1
%heun(A,B,fw,fw2,dt,t,x0)