Error using symengine dimensions do not match.

Hi, I am trying to use ode45 to solve a differential equation with the following codes (3 different codes)

First Function

function [A0,B,F,M ] = lagrange( )

%Orizume tis parametrus tu sistimatos me symbolics

syms I1 M1 R1 J2 J3 L4 I5 J8 M2 R2 J7 L10 K1 R1 K2 K3 C5 K8 K9 K7 R4 C9 ...

R10 V4 V6 T2 T3 G4 R4 T9 T8 T7 P5 P8 V10

%statheres parametroi systhmatos

syms V1 Q4 Q10 V1D V1DD TH THD THDD Q4D Q4DD V3D V5 V5D V5DD FD FDD Q10D Q10DD F

%opou:

%D=Proti paragwgos ton metavlitwn

%DD=Defteri paragwgos twn metavlitwn

syms PS VS %<--Exoterikes piges

%Va8moi Eleftherias

n=6;

%Anexartites metavlites defteris paragwgu

VariablesDD = [THDD FDD V1DD V5DD Q4DD Q10DD]

%Anexartites metavlites protis paragwgu

VariablesD = [THD FD V1D V5D Q4D Q10D]

%Anexartites metavlites

Variables = [TH F V1 V5 Q4 Q10 ]

%Kinitiki Energeia

T=0.5*(I1*(V1D)^2)+0.5*((M1/(R2)^2)*(THD)^2)+0.5*(J2*(THD)^2)+0.5*(J3*(THD)^2)+0.5*(L4*(Q4D)^2)+0.5*(I5*(V3D)^2)+0.5*(J8*(FD)^2)...

+0.5*((M2/(R2)^2)*(FD)^2)+0.5*(J7*(FD)^2)+0.5*(L10*(Q10D)^2);

%Dinamiki Energeia

U= ((1/2)*(K1/(R1^2))*TH^2)+((1/2)*K2*TH^2)+((1/2)*K3*TH^2)+((1/2)*(1/C5)*V5^2)+((1/2)*K8*F^2)+((1/2)*(K9/(R2^2))*F^2)+((1/2)*K7*F^2);

%Rayleigh

R=(0.5*R4*(Q4D^2))+(0.5*(C9/R2)*(FD^2))+(0.5*R10*(Q10D^2));

%Exoteriko Ergo

Wex= (PS*V1)-(V4*Q4)+(T2*TH)+(V6*Q4)-(T9*F)+(T3*TH)-(V5*P5)-(T8*F)+(V5*P8)+(T7*F)-(V10*Q10)-(VS*Q10);

%lagrange lisi.

L=T-U;

%dimiourgia pinakon pou tha xrisimopoiithoun stis eksisosis lagrange

eq=sym(zeros(n,1));

D=sym(zeros(n,1));

for i=1:n

D(i)=diff(Wex,Variables(i))-diff(R,VariablesD(i));

for j=1:n

eq(i)=eq(i)+diff(diff(L,VariablesD(i)),VariablesD(j))*VariablesDD(j);

end

eq(i)=eq(i)-diff(L,Variables(i))-D(i);

Fex(i,1)=diff((diff(Wex,PS)),Variables(i))*PS+...

diff(diff(Wex,VS),Variables(i))*VS;

end

%Anaparastasi ston xwro katastasis

[M]=equationsToMatrix([eq(1:n)], VariablesDD)

[C]=equationsToMatrix([eq(1:n)], VariablesD)

[K]=equationsToMatrix([eq(1:n)], Variables)

[ExoterikesPiges]=Fex

%Orizume tus voi8itikus pinakes(midenikos kai monadiaios)

O=sym(zeros(n)); I=sym(eye(n));

%Ipovivasmos Taxis

M(2,2) = 1; M(4,4)=1;

A0=[O I; -M^(-1)*K -M^(-1)*C]; A0=simplify(A0)

B=[O,O;O, M^(-1)]; B=simplify(B)

F=[0;0;0;0;0;0;0;0;0;Fex]; F=simplify(F)

disp('End Lagrange Method')

end

Second Function

function [ dydt ] = sol( t,x,A0,B,F )

syms PS VS %Prwta orizume ta symbolics

%Stin sinexeia antika8istume me times stis ekwterikes piges

F=subs(F,PS,20*cos(t));

F=subs(F,VS,40*cos(t));

%Stin sinexeia akolou8ei ipovivasmos taxis

dydt = A0*x+B*F;

dydt=double(dydt);

%dipli akrivia gia na sinadi me toys ipolipoys ipologismoys

end

Main Code

%Orizonte oi dinamikoi parametroi kai exoterikes piges

%sistimatos me symbolics

%dinamikes parametroi

syms I1 M1 R1 J2 J3 L4 I5 J8 M2 R2 J7 L10 K1 R1 K2 K3 C5 ...

K8 K9 K7 R4 C9 R10 V4 V6 T2 T3 G4 R4 T9 T8 T7 P5 P8 V10

syms VS PS %<--Exoterikes piges

%ARXIKOPIOISI

Dt=0.1; Tmax=100; t=0:Dt:Tmax; %Dt=Vima xronu ,Tmax=MAX XRONOS

InValues=[6;15;5;4;18;7;15;6;16;24;10;7]

%<--ARXIKES TIMES METAVLITON SISTIMATOS(12)

[A0,B,F] = lagrange(); %<--EPILISI LANGRANGE ME SIMBOLICS

%Epilisi sistimatos me tin methodo bond diagram

%run bond;

Elements={I1 M1 R1 J2 J3 L4 I5 J8 M2 R2 J7 L10 K1 R1 K2 K3 C5 K8 K9 K7 ...

R4 C9 R10 V4 V6 T2 T3 G4 R4 T9 T8 T7 P5 P8 V10 };

%(<--35)

%Dinume times stis pio panw parametrus

Values={ 15 30 10 20 30 3 10 15 25 27 10 15 15 10 25 35 26 28 20 5 0.5...

0.05 0.1 60 50 20 5 15 2 7 16 3 24 16 6};

%Antikatastasi ton statheron parametron me tis times toys

A0=subs(A0,Elements,Values);

B=subs(B,Elements,Values);

F=subs(F,Elements,Values);

%Epilisi diaforikwn exisosewn

disp('Epilisi diaforikon eksisoseon')

options = odeset('AbsTol',1e-4); %Orismos sfalmatos

[t,x]=ode45(sol,t,InValues,options,A0,B,F);

disp('telos epilisis diaforikon eksisoseon')

%Dimiurgia Grafikwn Parastasewn gia ka8e anexartiti metavliti se sxesi me

%ton xrono kai na emfanizontai oles se ena para8iro

subplot(2,3,1),plot(t,x(:,1)),title('TH vs t');

xlabel('Time (s)');ylabel('TH (rad)');

subplot(2,3,2);plot(t,x(:,2));title('F vs t');

xlabel('Time (s)');ylabel('F (rad)');

subplot(2,3,3);plot(t,x(:,3));title('V1 vs t')

xlabel('Time (s)');ylabel('V1 (m^3)')

subplot(2,3,4);plot(t,x(:,4));title('V5 vs t');

xlabel('Time (s)');ylabel('V5 (m^3)')

subplot(2,3,5);plot(t,x(:,5));title('Q4 vs t');

xlabel('Time (s)');ylabel('Q4 (C)')

subplot(2,3,6);plot(t,x(:,6));title('Q10 vs t')

xlabel('Time (s)');ylabel('Q10 (C)')

%Dimiurgia Grafikwn Parastasewn gia ka8e exoteriki pigi se sxesi me

%ton xrono kai na emfanizontai oles se ena para8iro

subplot(4,2,7),plot(t,140*cos(t)),title('PS vs t');

xlabel('Time (s)');ylabel('PS (Pa)');

subplot(4,2,8);plot(t,20*cos(t));title('VS vs t');

xlabel('Time (s)');ylabel('VS (Volts)');

and I am getting this error:

Error using symengine

Dimensions do not match.

Error in sym/privBinaryOp (line 1034)

Csym = mupadmex(op,args{1}.s, args{2}.s, varargin{:});

Error in * (line 329)

X = privBinaryOp(A, B, 'symobj::mtimes');

Error in sol (line 8)

dydt = A0*x+B*F;

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Does anyone knows how to solve this issue?