MATLAB code to use for the control of a harmonic oscillator.

function [x, u, Sf, tf] = ex3_3_5(a, b, r, x0) 
    % Control of a Harmonic Oscillator
    % Compute the solution to Riccati Equation
    [tb, S]=ode45(@fex3_3_5,[-10:0.1:0],zeros(3,1)); 
    % Compute Optimal Feedback Gains
    Sf=flipud(S); 
    tf=-flipud(tb); 
    K=-b/r*Sf(:,2:3); 
    x(:,1)=x0;
    u(1)=K(1,:)*x(:,1) ;
    % compute Closed-loop Response 
    for k=1:length(tf)-1,
        % Harmonic Oscillator System State Equations
        x(:,k+1) =expm((a+[0; b]*K(k,:))*(tf(k+1)-tf(k)))*x(:,k); 
        u(k+1)=K(k+1,:)*x(:,k+1);
    end
function sd=fex3_3_5(t,s)
    q=1*eye(2); om=0.8; del=0.1; b=1; r=1;
    sd =[-2*om^2*s(2)-b^2*s(2)^2+q(1, 1);
    s(1)-4*del*om*s(2)-om^2*s(3)-b^2*s(2)*s(3); 
    2*s(2)-4*del*om*s(3)-b^2*s(3)^2+q(2,2)];

please if anyone can help,

To output the following graph. Its from optimal control Frank Lewis.