Optimal control in Matlab

Hello, I need help with this Matlab code.

This is about an optimal control project.

It finds all the values until DV2 but it can't compute the last line: sol_h.

It gives no error but only the written "busy" and it gives no solutions (like it can't go on and stops).

May you tell me if you find any errors in initializing the variables or something else?

Thank you so much.

%EG1OC    Example 1 of optimal control tutorial.
%    This example is from D.E.Kirk's Optimal control theory: an introduction
%    example 5.1-1 on page 198 - 202

% State equations
%variables
syms U T R V1 V2 ;
%parameters
syms lambda k1 k2 k3 alpha1 alpha2 alpha3 gamma1 gamma2 phi1 phi2 delta1 delta2
%controls
syms u1 u2 u3 u4
dU = lambda - (1 + u1)*k1*U*V1 - (1 + u2)*k2*U*V2 - alpha1*U + gamma1*T + gamma2*R;
dT = (1 + u1)*k1*U*V1 - k3*T*V2 - alpha2*T - gamma1*T;
dR = (1 + u2)*k2*U*V2 - alpha3*R - gamma2*R + k3*T*V2;
dV1 = alpha2*T + gamma1*T + (1 + u3)*phi1*U + k3*T*V2 - delta1*V1;
dV2 = alpha3*R + gamma2*R + (1 + u4)*phi2*U - delta2*V2;

% Cost function inside the integral
syms g A1 A2 A3 A4 A5;
g = A1*U + (A2/2)*(u1)^2 + (A3/2)*(u2)^2 + (A4/2)*(u3)^2 + (A5/2)*(u4)^2;

% Hamiltonian
syms pU pT pR pV1 pV2 H;
H = g + pU*dU + pT*dT + pR*dR + pV1*dV1 + pV2*dV2;

% Costate equations
DpU = -diff(H,U);
DpT = -diff(H,T);
DpR = -diff(H,R);
DpV1 = -diff(H,V1);
DpV2 = -diff(H,V2);

% solve for control u
du1 = diff(H,u1);
sol_u1 = solve(du1, u1);
du2 = diff(H,u2);
sol_u2 = solve(du2, u2);
du3 = diff(H,u3);
sol_u3 = solve(du3, u3);
du4 = diff(H,u4);
sol_u4 = solve(du4, u4);

% Substitute u to state equations
DU = subs(dU, {u1 u2 u3 u4}, {sol_u1 sol_u2 sol_u3 sol_u4});
DT = subs(dT, {u1 u2 u3 u4}, {sol_u1 sol_u2 sol_u3 sol_u4});
DR = subs(dR, {u1 u2 u3 u4}, {sol_u1 sol_u2 sol_u3 sol_u4});
DV1 = subs(dV1, {u1 u2 u3 u4}, {sol_u1 sol_u2 sol_u3 sol_u4});
DV2 = subs(dV2, {u1 u2 u3 u4}, {sol_u1 sol_u2 sol_u3 sol_u4});

% convert symbolic objects to str ings for using 'dsolve'
eq1 = strcat('DU=',char(DU));
eq2 = strcat('DT=',char(DT));
eq3 = strcat('DR=',char(DR));
eq4 = strcat('DV1=',char(DV1));
eq5 = strcat('DV2=',char(DV2));
eq6 = strcat('DpU=',char(DpU));
eq7 = strcat('DpT=',char(DpT));
eq8 = strcat('DpR=',char(DpR));
eq9 = strcat('DpV1=',char(DpV1));
eq10 = strcat('DpV2=',char(DpV2));

sol_h = dsolve(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10);