Creating a system of function handle and manipulating it without cell array

I want to share a little bit of my knowledge because I had a rough time doing this simple thing. I wanted to solve a system of linear differential equation for a university homework using euler implicit method. The problem was that is didn't know how to :

1) Format the indidually defined differential function handle to create a system of equation.

2) Using the system of equation to get and format the non-linear system of equations required by the solver that our professor gave us (His solver is just a numerical method to find the root but work for systems of equations and vector of initial conditions)

So if you have the individually defined function handles like this for exemple :

t = 0 % Starting time
h = 0.01 % time step
y0 = [10 0] % Vector of Initial conditions for the system

% first order differential equations
f1 = @(t,y0) y0(2) ;
f2 = @(t,y0) -R*y0(2) - y0(1)/C  ;

1) To create a system of equation all you have to do is :

f = @(t,y0) [f1(t,y0) ; f2(t,y0) ]

so now you can evaluate a system of equation with a starting time t and a vector of initial conditions and it give your an output the same dimension of f

Euler implicit method requires you to to substitute the function with the unknown variables and , unlike Euler explicit method.

Euler explicit :

Euler implicit :

if we substitute by x

and we solve x for = 0 to find the roots with newtons or any other numerical method (this has to be done at each time step)

If we do that for every function we get the system of equation

2) As said earlier, I didn't know to to get from my system of first order differiential equations f to the non linear system of equation without having to type it all out by hand.

Turns out it was fairly simple. Our vector of initial conditions have to be the same length than our system f and so we use it to define the maximum value of i :

syms('x',[1 length(y0)]) % we define x to be substitued later

g= @(x)[] 
for i=1:1:length(y0)
k = @(x)[x(i)-y0(i)] 
g = @(x)[g(x); k(x)] %combining the functions
end
% at this point g has the same dimension as f and thus can be added and evaluated at f(t+h,x) to create the system of non linear equation
g = @(x) [g(x)-h*f(t+h,x)] % we assemble the functions to get the system g_system

So that's it, we have our system g and using the solver i get the root for each step.