Parametrised Implicit 1st order Euler Solver for ODEs

This project is an extension of the project shown here,

https://projects.skill-lync.com/projects/Solving-a-system-of-ODEs-using-Implicit-Euler-method-86502

In the above-shown project, I have used python to create an Implicit Euler ODE solver for 3 ODEs, If needed to solve for more or less than 3 ODEs, there needs to substantial changes that need to be made to the program, especially the Jacobian matrix,

This was the primary motivation to write the following code using MATLAB, using the same algorithm but to facilitate a more generalized solver,

In the following code, I have parametrized the Jacobian and the initial and final guess vectors of the Newton Rhapson Solver, to solve for "N" equations,

Hence, this prevents the tedious process of altering the Jacobian and the guess vectors for every case depending on the number of equations,

The code is given below,

clear all 
close all
clc


inty = [1,0,0];

[t,y1] = implicit_euler(inty,600,0.6);
t=transpose(t);

figure
hold on
plot(t,y1(:,1),'r')
plot(t,y1(:,2),'b')
plot(t,y1(:,3),'g')



function S = S1(y,ynew,h)

S(1) = y(1) + h*((-0.04*ynew(1)) + (10000*ynew(2)*ynew(3))) - ynew(1);
S(2) = y(2) + h*((0.04*ynew(1) - (10000*ynew(2)*ynew(3)) - (3*10000000*(ynew(2)^2)))) - ynew(2);
S(3) = y(3) + h*(3*10000000*(ynew(2)^2)) - ynew(3);

end

function J = numjacobian(y,ynew,h)

S = S1(y,ynew,h);
n = length(S);
diff = 1e-9;

for i = 1:n
    for j = 1:n
        
        d=zeros(n);
        d(j) = diff;
        
        for k = 1:n
            ynew1(k) = ynew(k) + d(k);
        end
        
        Jall = (S1(y,ynew1,h) - S1(y,ynew,h))/diff;
        J(i,j) = Jall(i);
        
        
    end
end

end

function Snewfinal = newtrhap(y,ynew,h)

S = S1(y,ynew,h); 

for i = 1:length(S)
    Sold(i) = ynew(i);
end


error = 9e9;
tol = 1e-9;
alpha = 1;

iter = 1;

while (error>tol)
    
    J = numjacobian(y,Sold,h);    
    F = S1(y,Sold,h);
    F = transpose(F);
    
    Snew = Sold - alpha*mtimes(inv(J),F);
    
    error = max(abs(Snew - Sold));
    
    Sold = Snew;
    
    Snewfinal = Snew(:,1);
    Snewfinal = transpose(Snewfinal);
    
    iter = iter + 1
    
end

end

function [t,Y] = implicit_euler(inty,tspan,dt)

t = 1:dt:tspan;
n = length(inty);

for ii = 1:n
    Y(1,ii) = inty(ii);
    yguess(ii) = 10;
end

for i = 2:length(t)
    
    Yall = newtrhap(Y((i-1),:),yguess,dt);
    Y(i,:) = Yall;
    yguess = Y(i,:);
    
end    
end

Here in this code, you need to define the input equations to the solver in the definition of Function S, here the ODEs are also supposed to be input in the Implicit Euler format as shown below,

Any "n" number of equations can be added to this solver it can solve for all the "n" variables without having to change anything else, the jacobian and guess vectors will automatically get replaced accordingly.

The result is also as consistent as the original solution using the hardcoded jacobian

SOLUTION USING PYTHON AND HARCODED JACOBIAN

SOLUTION USING MATLAB AND GENERALISED JACOBIAN