Solving a system of 2 ODEs with Implicit Euler Method? Please Help!

I have to find and plot the solution for this system of ODEs. Using ODE15s was easy, the hard part is that I must also solve this sytem using the implicit/backward euler method: dy1/dt = y(2); dy2/dt = 1000*(1-y(1)^2)*y(2)-y(1); t=0-3000 seconds, y1(0) = 2, y2(0) = 0.

The code for the euler method seems to work up until a certain time point (t = 807). After that, fsolve is not able to continue to find the right solution for some reason. Paste my code, run it and look at the plot. I have tried implementing other methods instead of fsolve, such as newtons method or successive substitution or even changing step size but i was unsuccessful in doing so. Please help me out!

dydt = @(t,y) [y(2);1000*(1-y(1)^2)*y(2) - y(1)]; %system of ODEs
t = 0:3000; %time vector
dt = t(2)-t(1); %time step
options = optimset('disp','off');

%Part A] Backward Euler
tic
yfe1 = 2;
yfe2 = 0;
yfe = [yfe1; yfe2];
yb(:,1) = yfe;
for i = 1:length(t)-1
    yfe(:,i+1) = yfe(:,i) + dydt(t(i),yfe(:,i))*dt; %explicit forward euler for initial guess
    fun = @(yn) yb(:,i) + dt*dydt(t(i+1),yn)-yn;  %set up function x-f(x)= g(x) = 0
    yb(:,i+1) = fsolve(fun,yb(:,i),options); %use fsolve to find next point
end
fprintf('\nThe Backward Euler method at t = 3000 gave y1,y2 = %f %f \n',yb(:,end))
toc

%Part B] ODE15 Method
tic
y0 = [2 0];
[tsolve, ysolve] = ode15s(@ode,t,y0);
tsolve = tsolve';
ysolve = ysolve';
fprintf('\nThe ODE15s method at t = 3000 gave y1,y2 = %f %f \n',ysolve(:,end))
toc

%plotting
hold on
plot(t,yb,'red')
plot(t,ysolve,'black')
xlabel('time [s]')
ylabel('y1, y2')
l1 = plot(nan, nan, 'r');
l2 = plot(nan, nan, 'black');
legend([l1, l2], 'Backward Euler', 'ODE15s')

function dy = ode(t,y)
dy(1) = y(2);
dy(2)= 1000*(1-y(1)^2)*y(2) - y(1);
dy = dy';
end

Answer 1

[Torsten]

As you can see from the plot obtained from ODE15S, there are discontinuities in the solution of the ODE.

I don't think you will be able to catch these points with your Euler backwards.

But the equations you try to solve with "fsolve" can be solved analytically. Although it is quite of cheating (because in the usual cases, you won't be able to solve the nonlinear equations analytically), inserting this analytical solution for yb(:,i+1) should help.

Use symbolic math to calculate the analytical solution for yb(:,i+1) once and copy the formula thus obtained into a function file. This function can subsequently be called from your program.