% matlab script to compare FE, BE, trapezoidal rule on lambdaSin
% ODE

clear all, close all

lambda = -100;
t0 = 0;
dt = 0.01;
y0 = 1;
tf = 3;

t = [t0:dt:tf];
nSteps = length(t)-1;
yFE = [y0; zeros(nSteps,1)];
yBE = [y0; zeros(nSteps,1)];
yTR = [y0; zeros(nSteps,1)];

% FE
for n=1:nSteps-1,
  yFE(n+1) = yFE(n) + dt*lambda*(yFE(n)-sin(t(n)));
end

figure(1)
plot(t,yFE)

% BE
for n=1:nSteps-1,
  yBE(n+1) = (yBE(n) + dt*lambda*(-sin(t(n+1))))/(1-dt*lambda);
end

figure(2)
plot(t,yBE)

% TR
for n=1:nSteps-1,
  yTR(n+1) = (yTR(n) + dt*lambda/2*(yTR(n)-sin(t(n))-sin(t(n+1))))/(1-dt*lambda/2);
end
figure(3)
plot(t,yTR)