issue in solving set of ode's

Here i gave my script with dummy R, G, L, g matrices(because of confidentiality).

Here iam forming the matrices G, L and g interms of theta and extracting these variable matrices to solve set of odes. In that iam substituting 26th solution(y(26)) in place of theta. I attached photo better look into that.

I want to know whether the steps from "commonvars" are correct or not.

syms phi theta

a = 2*phi+5;
b = diff(a);

R = rand(27, 27);

Gi = b*cos(theta)*rand(25, 25);

G = [Gi zeros(25, 2);
    zeros(2, 27)]

Li = sin(theta+5)*rand(25, 25);
 
L = [Li zeros(25, 2);
    zeros(2, 25) [1 0; 0 1]]

g = (sin(theta)+cos(theta+2))*rand(25, 25)

commonvars = unique([symvar(G), symvar(L), symvar(g)]);  %probably just theta
Gfun = matlabFunction(G, 'vars', commonvars);
Lfun = matlabFunction(L, 'vars', commonvars);
gfun = matlabFunction(g, 'vars', commonvars);

tspan = [0 1];
myfun = @(t,y)scriptname(t, y, R, Gfun, Lfun, gfun);
y0 = zeros(27, 1);
[t, y] = ode45(myfun, tspan, y0);

h = figure;
% plot 
for i=1:25
 plot(t,y(:, i));
end
   
function dydt = scriptname(t, y, R, Gfun, Lfun, gfun)
    f = 60;  % frequency
    Wr = 2*pi*f;  % Angular velocity of rotor
    inertia = 0.05;
    % evaluation of g (numerical) with theta = y(27)
    gn = gfun(y(27))
    for i=1:25
        I(i,1)=y(i);
    end
    Te = 1/2*I'*gn*I
    
    V=[sqrt(2)*400/sqrt(3)*cos(Wr*t);
    sqrt(2)*400/sqrt(3)*cos(Wr*t+2.*pi/3.);
    sqrt(2)*400/sqrt(3)*cos(Wr*t-2.*pi/3.);
    zeros(Nr+1, 1);
    0;
    Te/inertia;
    y(26)];
    
    % evaluation of G and L (numerical) with theta = y(27)
    Gn = Gfun(y(27))
    Ln = Lfun(y(27))
    dydt = Ln\(V-R*y-y(26)*Gn*y)
end

Actually the code was running fine and giving waveforms also but not as expected.

I want to make sure whether steps for ode calculation are correct or not.