Rozwiązanie układu równań różniczkowych dla transmitancji

Mam za zadanie napisać program który rozwiąże układ równań różniczkowych opisany transmitancją:

G(s)=k/(s*(Tws+1)(Tms+1)+k)

Dostałem za przykład taki kod:

#include <iostream.h>
#include <conio.h>
#include <stdio.h>
#include <math.h>
/*program rozwiazuje uklad dwoch rownan rozniczkowych powstaly z rownania
rozniczkowego drugiego rzedu czlonu oscylacyjnego metoda trapezow
	 T^2*d^2y/dy^2+2*ksi*T*dy/dt+y=k*f(t)
przez podstawienie y=y1,dy1/dt=y2. W wyniku otrzymamy uklad rownan pierwszego
 rzedu
    dy1/dy=y2
    dy2/dt=-2*ksi/T*y2-1/T/T*y1+k/T/T*f(t)

Oznaczenia k -wspolczynnik wzmocnienia;
	   ksi - wspolczynnik tlumienia
	   T - stala czasowa
nalezy podstawic warunki poczatkowe y1(t=t0) i y2(t=t0)
*/
const float pi=M_PI;
//---------------------------------
 float f(float t)
  {
    float f0=100.0; //czestotliwosc
   return 1.0;
    //return sin(2*f0*pi*t);
  }

 int main()
  {
   FILE *wsk; //wskaznik typu FILE
   wsk=fopen("d:\\1\\rozw1.dat","w");
   float y1,y2,dy1,dy2,t,dt,tk,g,ksi,T,k;
   int n,m;
   n=2; m=2;
   float**ptrArrayY; //macierz wskaznikow calkowanych wartosi
                       //wskaznikow zmiennych calkowanych Y)
   float**ptrArraydY;  //tablica  wskaznikow pochodnych)
   float* ArraydY_1;// tablica pierwszych wartosci dy1 i dy2;
   //Dynamiczny przydzial pamieci
  if( (ptrArrayY  = new float*[n])==NULL ||
      (ptrArraydY = new float * [n])==NULL ||
      (ArraydY_1 =  new float [n])==NULL )
  {
    printf(" ERROR: Brak pamieci dla macierzy");
    cout<<endl;
    exit(1);
  }
   t=0.0;  dt=0.001;  tk=2.0; k=1.0; ksi=0.3; T=0.1;
   y1=0.0; y2=0.0;
   y1  = 0.0;   ptrArrayY[0] = & y1;
   g=f(t);
   dy1=y2;  ptrArraydY[0] = &dy1;
   y2   = 0.0;   ptrArrayY[1] = & y2;
   dy2=-2*ksi/T*y2-1/T/T*y1+k/T/T*f(t);  ptrArraydY[1] = &dy2;
  // *(ArraydY_1)=0.0;

      for (int i=0; i < n; i++)
       *(ArraydY_1+i)=*ptrArraydY[i];

     do
     {
      t+=dt;
      g=f(t);
      dy1=y2;
      dy2=-2*ksi/T*y2-1/T/T*y1+k/T/T*f(t);
       for (int i=0; i<n; i++)
    *(ArraydY_1+i)=*ptrArraydY[i];
      for( int i=0; i < n; i++ )
     *ptrArrayY[i] += (*ptrArraydY[i] + *(ArraydY_1+i)) * dt/2;
      fprintf(wsk,"\n %f %f %f",t,*ptrArrayY[0],*ptrArrayY[1]);
     // *(ArraydY_1+i))
     } while(t<tk);

      fclose(wsk);
      getch();
      return 0;
  }

Po wielu próbach odtworzenia tego programu dla równania o wyższym rzędzie, przychodzę prosić o jakieś ukierunkowanie. To do czego sam doszedłem to układ równań:
link

Nie bardzo wiem jak uzależnić kolejne pochodne tak aby zostały same równania pierwszego rzędu.