Optimal synthesis control problem

I want to solve a synthesis problem in control theory:

By Pontryagin's maximum principle I've getten:

The dimensions of matrix are big, for example

I want to write a program in MATLAB that will find me the optimal control and dependence of variables on time.

The algorithm is as follows: take a random variable:

after take a time from for which . Is my reasoning correct or are any numerical methods which can solve the synthesizing problem?

I tried to write the MATLAB code, but I'm sure It's wrong.

Thanks for the help!

clc;
clear;

Ax = 0.25; %kg/s
Ay = 0.25; %kg/s
Az = 0.25; %kg/s
b = 1.140*10^(-7);
g = 9.8; % m/s^2;
k = 2.980*10^(-7);
l = 0.225;
u_c = 1;
m = 0.468; %kgs
Ixx = 4.856*10^(-3);
Iyy = 4.856*10^(-3);
Izz = 8.801*10^(-3);

A = [0   1   0   0   0   0    0   0   0   0   0   0;
     0 -Ax/m 0   0   0   0    0   0   g   0   0   0;
     0   0   0   1   0   0    0   0   0   0   0   0;
     0   0   0 -Ay/m 0   0   -g   0   0   0   0   0;
     0   0   0   0   0   1    0   0   0   0   0   0;
     0   0   0   0   0 -Az/m  0   0   0   0   0   0;
     0   0   0   0   0   0    0   1   0   0   0   0;
     0   0   0   0   0   0    0   0   0   0   0   0;
     0   0   0   0   0   0    0   0   1   0   0   0;
     0   0   0   0   0   0    0   0   0   0   0   0;
     0   0   0   0   0   0    0   0   0   0   0   1;
     0   0   0   0   0   0    0   0   0   0   0   0;
     ];
B = [ 0     0       0      0;
      0     0       0      0;
      0     0       0      0;
      0     0       0      0;
      0     0       0      0;
     k/m   k/m     k/m    k/m;
      0     0       0      0;
      0  -k*l/Ixx   0   k*l/Ixx;
      0     0       0      0;
   -k*l/Iyy 0    -k*l/Iyy  0;
       0     0      0      0;
    b/Izz -b/Izz  b/Izz -b/Izz;
      ];

%initial state Psi_t0 random
Psi_t0 = [0.03   0.1   0.4   0   0   0.3    0   0.1   0   0.1   0   1].';

tspan = [0:1:100];

[t, Psi] = ode45(@(t, Psi) -A.'*Psi, tspan, Psi_t0);

Psi = Psi.';
%control
u1 = u_c*sign(B.'*Psi);

x_0 = [0.03   0.1   0.4   0   0   0.3    0   0.1   0   0.1   0   1].';
for i=1:1:100
    tspan =[i-1 i];
    u = u1(:,i);
    [t, x] = ode45(@(t, x) A*x + B*u, tspan, x_0);
    x_0 = x(end, :);
end