I think the ode45 code that I wrote works, but how do I plot it together with both forward and backward euler methods?

% values of a & b
a = 7;
b = 3;
L = 0.01*b;
R = 10*a;
C = (0.0001)/a;
t = [0;((40*L)/R)];
f = 1/t;
V = 0;

omega = 2*pi*f;

%%% First-order equations
%%% dI = J;
%%% d2I =dJ;
%%% dJ =((V1*omega)/L)*cos(omega*t)-(R/L)*J-(I/(L*C));
%%% Note: R,C & L are considered constants

% Boundary Conditions
ua = 0;
ub = 0;

% Solve system with initital conditions through ode45() solver
mu = 1;
[T,Y] = ode45(@(t,v) [v(2);mu*(1-v(1).^2).*v(2)-v(1)],[0 ((40*L)/R)],[ua;ub]);

% h value given
h = 0.001;
N = round(((40*L)/R)/h);

%%% initial conditions for Euler methods are the same
%%% as boundary conditions

%initial conditions for Euler Method
xf(1) = ua;
yf(1) = ub;
tf(1) = 0;

% Forward Euler method
for n = 1:N
    tf(n+1)=tf(n)+h;
    xf(n+1)=xf(n)+yf(n)*h;
    yf(n+1)=yf(n)+(mu*(1-xf(n)^2)*yf(n)-xf(n))*h;
end

% Backward Euler method
xb(1) = ua;
yb(1) = ub;
tb(1) = 0;

% Backward Euler method
for n = 1:N
  tb(n+1) = tb(n)+h;
  fun = @(u,v)[(u-xb(n))/h - v;(v-yb(n))/h - (mu*(1-u^2)*v-u)];
  sol = fsolve(@(x)fun(x(1),x(2)),[xb(n),yb(n)],optimset('Display','none'));
  xb(n+1) = sol(1);
  yb(n+1) = sol(2);
end
% Plotting
figure
hold on
plot(T,Y(:,1))
plot(tf,xf)
plot(tb,xb)
hold off

Answer 1

[Torsten]

I changed your code as to make it work (see above).