clear vb vf vt t;
lambda = -1000;
Tmax   = 10;
dt = Tmax/N;

t = linspace(0,Tmax,N+1);
vf = zeros(1,N+1); % It is important to allocate the space for
                   % speed as N gets large.
vb = zeros(1,N+1);
vt = zeros(1,N+1);

% Exact solution
omega = 1;
c2 = 100/(omega^2 + lambda^2);
c1 = 1 + omega*c2;
ue = c1*exp(lambda*t) + c2*(-lambda*sin(omega*t) - omega*cos(omega*t));


% Forward Euler method
vf(1) = 1.0;
for n = 1:N,
  vf(n+1) = vf(n) + dt*(lambda*vf(n) + 100*sin(t(n)));
end

% Backward Euler method
vb(1) = 1.0;
for n = 1:N,
  vb(n+1) = (vb(n) + dt*100*sin(t(n+1)))/(1 - dt*lambda);
end

% Trapezoidal method
vt(1) = 1.0;
for n = 1:N,
  vt(n+1) = (vt(n) + 0.5*dt*lambda*vt(n) + 0.5*dt*100*(sin(t(n))+sin(t(n+1))))/(1 - 0.5*dt*lambda);
end

%plot(t,ue,'r');
%hold on;
%plot(t,vt);
%axis([0,10,-1,1]);
%hold off;