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

Can any one help with matlab code for a harmonic osillator

function [x, u, Sf, tf] = ex3_3_5(a, b, r, x0)

  om=0.8; del=0.1; b=1; r=1;
  A= [0 1;-om^2 -2*del*om]; B= [0;b];
  C= [1 0]; D= 0;
  Q= [0.1 0;0 0.1]; R= 1;

  % Control of a Harmonic Oscillator
  % Compute the solution to Riccati Equation
  [tb, S]=ode45(@HW3_4,[-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

  subplot 221
  plot(t,y(:,1),'b')
  xlabel('time[s]');
  legend('q=0.01','q=0.1','q=1');
  title('state x_1');
  hold on
   
  subplot 222
  plot(t,u,'b')
  xlabel('time[s]');
  legend('q=0.01','q=0.1','q=1');
  title('control u');
  hold on
  
  subplot 223
  range=[0 10];
  ic=[0;0;0];
  [t, S]=ode45(@HW3_4,range,ic);
  plot(t,S(:,1),'k-',t,S(:,2),'b',t,S(:,3),'r')
  legend('s(1)','s(2)','s(3)')
  xlabel('time[s]');
  title('q=0.01');
  hold on
  
  subplot 224
  range=[0 10];
  ic=[0;0;0];
  [t, S]=ode45(@HW3_4,range,ic);
  plot(t,S(:,1),'k-',t,S(:,2),'b',t,S(:,3),'r')
  legend('s(1)','s(2)','s(3)')
  xlabel('time[s]');
  title('q=1');

function sd=HW3_4(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)];