PHI = @(t) [ 1, t, t/3-(2*exp(-3*t))/9+2/9, (2*t)/3+(2*exp(-3*t))/9-2/9;
0, 1, (5*exp(-3*t))/12-(3*exp(t))/4+1/3, 2/3-exp(t)/4-(5*exp(-3*t))/12;
0, 0, exp(-3*t)/4+(3*exp(t))/4, exp(t)/4-exp(-3*t)/4;
0, 0, (3*exp(t))/4-(3*exp(-3*t))/4, (3*exp(-3*t))/4+exp(t)/4];
PHIT = @(t) transpose(PHI(t));
B=[0;1;2;1];
BT=transpose(B);
GRAMi = @(t) PHI(-t)*B*BT*PHIT(-t);
GRAMfinal = @(t) integral(GRAMi, 0, t, 'ArrayValued',1)
%InGRAMfinal = @(t) inv(GRAMfinal(t));
X0=[0;2;0;1];
t=1;
A= GRAMfinal(t)
B=inv(A)
UT= @(t) -(BT*PHIT(-t)*B(1)*X0);
TUT = @(t) 2*(-(BT*PHIT(-t)*B(1)*X0));
syms x1(t) x2(t) x3(t) x4(t)
ode1 = diff(x1) == 1*x2 + 1*x3;
ode2 = diff(x2) == -2*x3 + 1*x4 +UT;
ode3 = diff(x3) == 1*x4 + TUT;
ode4 = diff(x4) == 3*x3 + (-2)*x4 + UT;
odes = [ode1; ode2; ode3; ode4]
S = dsolve(odes)
x1Sol(t) = S.x1
x2Sol(t) = S.x2
x3Sol(t) = S.x3
x4Sol(t) = S.x4
cond1 = x1(0) == 0;
cond2 = x2(0) == 2;
cond3 = x3(0) == 0;
cond4 = x4(0) == 1;
conds = [cond1; cond2; cond3; cond4];
[x1Sol(t), x2Sol(t),x3Sol(t) , x4Sol(t)] = dsolve(odes,conds)
t=0:0.001:1;
grid on
hold on
plot(t,x1Sol(t),'black', 'linewidth' ,2)
hold on
plot(t,x2Sol(t), 'red','linewidth' ,2)
hold on
plot(t,x3Sol(t), 'black','linewidth' ,2)
hold on
plot(t,x4Sol(t),'red','linewidth' ,2)
hold on
xlabel('t [0,1]')
ylabel('Four components of the state x1, x2, x3 and x4')
