function g=makeTs(t,T)

g1=exp(T)./(-1 + exp(T));
g2=(exp(T).*((-1 + exp(T)).*t + T))./(-1 + exp(T)).^2;
g3=(exp(T).*((-1 + exp(T)).^2.*t.^2 + 2.*(-1 + exp(T)).*t.*T + (1 + exp(T)).*T.^2))./(-1 + exp(T)).^3;
g4=(exp(T).*((-1 + exp(T)).^3.*t.^3 + 3.*(-1 + exp(T)).^2.*t.^2.*T + 3.*(-1 + exp(2.*T)).*t.*T.^2 + (1 + 4.*exp(T) + exp(2.*T)).*T.^3))./(-1 + exp(T)).^4;
g5=(exp(T).*((-1 + exp(T)).^4.*t.^4 + 4.*(-1 + exp(T)).^3.*t.^3.*T + 6.*(-1 + exp(T)).^2.*(1 + exp(T)).*t.^2.*T.^2 + 4.*(-1 - 3.*exp(T) + 3.*exp(2.*T) + exp(3.*T)).*t.*T.^3 + ...
   (1 + 11.*exp(T) + 11.*exp(2.*T) + exp(3.*T)).*T.^4))./(-1 + exp(T)).^5;
g6=1./(-1 + exp(T)).^6.*exp(T).*((-1 + exp(T)).^5.*t.^5 + 5.*(-1 + exp(T)).^4.*t.^4.*T + 10.*(-1 + exp(T)).^3.*(1 + exp(T)).*t.^3.*T.^2 + ...
    10.*(-1 + exp(T)).^2.*(1 + 4.*exp(T) + exp(2.*T)).*t.^2.*T.^3 + 5.*(-1 - 10.*exp(T) + 10.*exp(3.*T) + exp(4.*T)).*t.*T.^4 + ...
    (1 + 26.*exp(T) + 66.*exp(2.*T) + 26.*exp(3.*T) + exp(4.*T)).*T.^5);
g7=1./(-1 + exp(T)).^7.*exp(T).*((-1 + exp(T)).^6.*t.^6 + 6.*(-1 + exp(T)).^5.*t.^5.*T + 15.*(-1 + exp(T)).^4.*(1 + exp(T)).*t.^4.*T.^2 + ...
    20.*(-1 + exp(T)).^3.*(1 + 4.*exp(T) + exp(2.*T)).*t.^3.*T.^3 + 15.*(-1 + exp(T)).^2.*(1 + 11.*exp(T) + 11.*exp(2.*T) + exp(3.*T)).*t.^2.*T.^4 + ... 
    6.*(-1 - 25.*exp(T) - 40.*exp(2.*T) + 40.*exp(3.*T) + 25.*exp(4.*T) + exp(5.*T)).*t.*T.^5 + ...
    (1 + 57.*exp(T) + 302.*exp(2.*T) + 302.*exp(3.*T) + 57.*exp(4.*T) + exp(5.*T)).*T.^6);
g8=1./(-1 + exp(T)).^8.*exp(T).*((-1 + exp(T)).^7.*t.^7 + 7.*(-1 + exp(T)).^6.*t.^6.*T + 21.*(-1 + exp(T)).^5.*(1 + exp(T)).*t.^5.*T.^2 + ...
    35.*(-1 + exp(T)).^4.*(1 + 4.*exp(T) + exp(2.*T)).*t.^4.*T.^3 + 35.*(-1 + exp(T)).^3.*(1 + 11.*exp(T) + 11.*exp(2.*T) + exp(3.*T)).*t.^3.*T.^4 + ... 
    21.*(-1 + exp(T)).^2.*(1 + 26.*exp(T) + 66.*exp(2.*T) + 26.*exp(3.*T) + exp(4.*T)).*t.^2.*T.^5 + ...
    7.*(-1 - 56.*exp(T) - 245.*exp(2.*T) + 245.*exp(4.*T) + 56.*exp(5.*T) + exp(6.*T)).*t.*T.^6 + ...
    (1 + 120.*exp(T) + 1191.*exp(2.*T) + 2416.*exp(3.*T) + 1191.*exp(4.*T) + 120.*exp(5.*T) + exp(6.*T)).*T.^7);
g9=1./(-1 + exp(T)).^9.*exp(T).*((-1 + exp(T)).^8.*t.^8 + 8.*(-1 + exp(T)).^7.*t.^7.*T + 28.*(-1 + exp(T)).^6.*(1 + exp(T)).*t.^6.*T.^2 + ...
    56.*(-1 + exp(T)).^5.*(1 + 4.*exp(T) + exp(2.*T)).*t.^5.*T.^3 + 70.*(-1 + exp(T)).^4.*(1 + 11.*exp(T) + 11.*exp(2.*T) + exp(3.*T)).*t.^4.*T.^4 + ...
    56.*(-1 + exp(T)).^3.*(1 + 26.*exp(T) + 66.*exp(2.*T) + 26.*exp(3.*T) + exp(4.*T)).*t.^3.*T.^5 + ...
    28.*(-1 + exp(T)).^2.*(1 + 57.*exp(T) + 302.*exp(2.*T) + 302.*exp(3.*T) + 57.*exp(4.*T) + exp(5.*T)).*t.^2.*T.^6 + ...
    8.*(-1 - 119.*exp(T) - 1071.*exp(2.*T) - 1225.*exp(3.*T) + 1225.*exp(4.*T) + 1071.*exp(5.*T) + 119.*exp(6.*T) + exp(7.*T)).*t.*T.^7 + ...
    (1 + 247.*exp(T) + 4293.*exp(2.*T) + 15619.*exp(3.*T) + 15619.*exp(4.*T) + 4293.*exp(5.*T) + 247.*exp(6.*T) + exp(7.*T)).*T.^8);
g10=1./(-1 + exp(T)).^10.*exp(T).*((-1 + exp(T)).^9.*t.^9 + 9.*(-1 + exp(T)).^8.*t.^8.*T + 36.*(-1 + exp(T)).^7.*(1 + exp(T)).*t.^7.*T.^2 + ...
    84.*(-1 + exp(T)).^6.*(1 + 4.*exp(T) + exp(2.*T)).*t.^6.*T.^3 + 126.*(-1 + exp(T)).^5.*(1 + 11.*exp(T) + 11.*exp(2.*T) + exp(3.*T)).*t.^5.*T.^4 + ...
    126.*(-1 + exp(T)).^4.*(1 + 26.*exp(T) + 66.*exp(2.*T) + 26.*exp(3.*T) + exp(4.*T)).*t.^4.*T.^5 + ...
    84.*(-1 + exp(T)).^3.*(1 + 57.*exp(T) + 302.*exp(2.*T) + 302.*exp(3.*T) + 57.*exp(4.*T) + exp(5.*T)).*t.^3.*T.^6 + ...
    36.*(-1 + exp(T)).^2.*(1 + 120.*exp(T) + 1191.*exp(2.*T) + 2416.*exp(3.*T) + 1191.*exp(4.*T) + 120.*exp(5.*T) + exp(6.*T)).*t.^2.*T.^7 + ...
    9.*(-1 + exp(T)).*(1 + 247.*exp(T) + 4293.*exp(2.*T) + 15619.*exp(3.*T) + 15619.*exp(4.*T) + 4293.*exp(5.*T) + 247.*exp(6.*T) + exp(7.*T)).*t.*T.^8 + ...
    (1 + 502.*exp(T) + 14608.*exp(2.*T) + 88234.*exp(3.*T) + 156190.*exp(4.*T) + 88234.*exp(5.*T) + 14608.*exp(6.*T) + 502.*exp(7.*T) + exp(8.*T)).*T.^9);
g11=1./(-1 + exp(T)).^11.*exp(T).*((-1 + exp(T)).^10.*t.^10 + 10.*(-1 + exp(T)).^9.*t.^9.*T + 45.*(-1 + exp(T)).^8.*(1 + exp(T)).*t.^8.*T.^2 + ...
    120.*(-1 + exp(T)).^7.*(1 + 4.*exp(T) + exp(2.*T)).*t.^7.*T.^3 + 210.*(-1 + exp(T)).^6.*(1 + 11.*exp(T) + 11.*exp(2.*T) + exp(3.*T)).*t.^6.*T.^4 + ...
    252.*(-1 + exp(T)).^5.*(1 + 26.*exp(T) + 66.*exp(2.*T) + 26.*exp(3.*T) + exp(4.*T)).*t.^5.*T.^5 + ...
    210.*(-1 + exp(T)).^4.*(1 + 57.*exp(T) + 302.*exp(2.*T) + 302.*exp(3.*T) + 57.*exp(4.*T) + exp(5.*T)).*t.^4.*T.^6 + ...
    120.*(-1 + exp(T)).^3.*(1 + 120.*exp(T) + 1191.*exp(2.*T) + 2416.*exp(3.*T) + 1191.*exp(4.*T) + 120.*exp(5.*T) + exp(6.*T)).*t.^3.*T.^7 + ...
    45.*(-1 + exp(T)).^2.*(1 + 247.*exp(T) + 4293.*exp(2.*T) + 15619.*exp(3.*T) + 15619.*exp(4.*T) + 4293.*exp(5.*T) + 247.*exp(6.*T) + exp(7.*T)).*t.^2.*T.^8 + ...
    10.*(-1 + exp(T)).*(1 + 502.*exp(T) + 14608.*exp(2.*T) + 88234.*exp(3.*T) + 156190.*exp(4.*T) + 88234.*exp(5.*T) + 14608.*exp(6.*T) + 502.*exp(7.*T) + ...
      exp(8.*T)).*t.*T.^9 + (1 + 1013.*exp(T) + 47840.*exp(2.*T) + 455192.*exp(3.*T) + 1310354.*exp(4.*T) + 1310354.*exp(5.*T) + 455192.*exp(6.*T) + ...
      47840.*exp(7.*T) + 1013.*exp(8.*T) + exp(9.*T)).*T.^10);
g12=1./(-1 + exp(T)).^12.*exp(T).*((-1 + exp(T)).^11.*t.^11 + 11.*(-1 + exp(T)).^10.*t.^10.*T + 55.*(-1 + exp(T)).^9.*(1 + exp(T)).*t.^9.*T.^2 + ...
    165.*(-1 + exp(T)).^8.*(1 + 4.*exp(T) + exp(2.*T)).*t.^8.*T.^3 + 330.*(-1 + exp(T)).^7.*(1 + 11.*exp(T) + 11.*exp(2.*T) + exp(3.*T)).*t.^7.*T.^4 + ...
    462.*(-1 + exp(T)).^6.*(1 + 26.*exp(T) + 66.*exp(2.*T) + 26.*exp(3.*T) + exp(4.*T)).*t.^6.*T.^5 + ...
    462.*(-1 + exp(T)).^5.*(1 + 57.*exp(T) + 302.*exp(2.*T) + 302.*exp(3.*T) + 57.*exp(4.*T) + exp(5.*T)).*t.^5.*T.^6 + ...
    330.*(-1 + exp(T)).^4.*(1 + 120.*exp(T) + 1191.*exp(2.*T) + 2416.*exp(3.*T) + 1191.*exp(4.*T) + 120.*exp(5.*T) + exp(6.*T)).*t.^4.*T.^7 + ...
    165.*(-1 + exp(T)).^3.*(1 + 247.*exp(T) + 4293.*exp(2.*T) + 15619.*exp(3.*T) + 15619.*exp(4.*T) + 4293.*exp(5.*T) + 247.*exp(6.*T) + exp(7.*T)).*t.^3.*T.^8 + ...
    55.*(-1 + exp(T)).^2.*(1 + 502.*exp(T) + 14608.*exp(2.*T) + 88234.*exp(3.*T) + 156190.*exp(4.*T) + 88234.*exp(5.*T) + 14608.*exp(6.*T) + 502.*exp(7.*T) + ...
      exp(8.*T)).*t.^2.*T.^9 + 11.*(-1 + exp(T)).*(1 + 1013.*exp(T) + 47840.*exp(2.*T) + 455192.*exp(3.*T) + 1310354.*exp(4.*T) + 1310354.*exp(5.*T) + ...
      455192.*exp(6.*T) + 47840.*exp(7.*T) + 1013.*exp(8.*T) + exp(9.*T)).*t.*T.^10 + ...
    (1 + 2036.*exp(T) + 152637.*exp(2.*T) + 2203488.*exp(3.*T) + 9738114.*exp(4.*T) + 15724248.*exp(5.*T) + 9738114.*exp(6.*T) + 2203488.*exp(7.*T) + ...
    152637.*exp(8.*T) + 2036.*exp(9.*T) + exp(10.*T)).*T.^11);
 g13=1./(-1 + exp(T)).^13.*exp(T).*((-1 + exp(T)).^12.*t.^12 + 12.*(-1 + exp(T)).^11.*t.^11.*T + 66.*(-1 + exp(T)).^10.*(1 + exp(T)).*t.^10.*T.^2 + ...
    220.*(-1 + exp(T)).^9.*(1 + 4.*exp(T) + exp(2.*T)).*t.^9.*T.^3 + 495.*(-1 + exp(T)).^8.*(1 + 11.*exp(T) + 11.*exp(2.*T) + exp(3.*T)).*t.^8.*T.^4 + ...
    792.*(-1 + exp(T)).^7.*(1 + 26.*exp(T) + 66.*exp(2.*T) + 26.*exp(3.*T) + exp(4.*T)).*t.^7.*T.^5 + ...
    924.*(-1 + exp(T)).^6.*(1 + 57.*exp(T) + 302.*exp(2.*T) + 302.*exp(3.*T) + 57.*exp(4.*T) + exp(5.*T)).*t.^6.*T.^6 + ...
       792.*(-1 + exp(T)).^5.*(1 + 120.*exp(T) + 1191.*exp(2.*T) + 2416.*exp(3.*T) + 1191.*exp(4.*T) + 120.*exp(5.*T) + exp(6.*T)).*t.^5.*T.^7 + ...
    495.*(-1 + exp(T)).^4.*(1 + 247.*exp(T) + 4293.*exp(2.*T) + 15619.*exp(3.*T) + 15619.*exp(4.*T) + 4293.*exp(5.*T) + 247.*exp(6.*T) + exp(7.*T)).*t.^4.*T.^8 + ...
    220.*(-1 + exp(T)).^3.*(1 + 502.*exp(T) + 14608.*exp(2.*T) + 88234.*exp(3.*T) + 156190.*exp(4.*T) + 88234.*exp(5.*T) + 14608.*exp(6.*T) + 502.*exp(7.*T) + ...
      exp(8.*T)).*t.^3.*T.^9 + 66.*(-1 + exp(T)).^2.*(1 + 1013.*exp(T) + 47840.*exp(2.*T) + 455192.*exp(3.*T) + 1310354.*exp(4.*T) + 1310354.*exp(5.*T) + ...
      455192.*exp(6.*T) + 47840.*exp(7.*T) + 1013.*exp(8.*T) + exp(9.*T)).*t.^2.*T.^10 + ...
    12.*(-1 + exp(T)).*(1 + 2036.*exp(T) + 152637.*exp(2.*T) + 2203488.*exp(3.*T) + 9738114.*exp(4.*T) + 15724248.*exp(5.*T) + 9738114.*exp(6.*T) + ...
      2203488.*exp(7.*T) + 152637.*exp(8.*T) + 2036.*exp(9.*T) + exp(10.*T)).*t.*T.^11 + ...
    (1 + 4083.*exp(T) + 478271.*exp(2.*T) + 10187685.*exp(3.*T) + 66318474.*exp(4.*T) + 162512286.*exp(5.*T) + 162512286.*exp(6.*T) + ...
    66318474.*exp(7.*T) + 10187685.*exp(8.*T) + 478271.*exp(9.*T) + 4083.*exp(10.*T) + exp(11.*T)).*T.^12);
 g14=1./(-1 + exp(T)).^14.*exp(T).*((-1 + exp(T)).^13.*t.^13 + 13.*(-1 + exp(T)).^12.*t.^12.*T + 78.*(-1 + exp(T)).^11.*(1 + exp(T)).*t.^11.*T.^2 + ...
    286.*(-1 + exp(T)).^10.*(1 + 4.*exp(T) + exp(2.*T)).*t.^10.*T.^3 + 715.*(-1 + exp(T)).^9.*(1 + 11.*exp(T) + 11.*exp(2.*T) + exp(3.*T)).*t.^9.*T.^4 + ...
    1287.*(-1 + exp(T)).^8.*(1 + 26.*exp(T) + 66.*exp(2.*T) + 26.*exp(3.*T) + exp(4.*T)).*t.^8.*T.^5 + ...
    1716.*(-1 + exp(T)).^7.*(1 + 57.*exp(T) + 302.*exp(2.*T) + 302.*exp(3.*T) + 57.*exp(4.*T) + exp(5.*T)).*t.^7.*T.^6 + ...
    1716.*(-1 + exp(T)).^6.*(1 + 120.*exp(T) + 1191.*exp(2.*T) + 2416.*exp(3.*T) + 1191.*exp(4.*T) + 120.*exp(5.*T) + exp(6.*T)).*t.^6.*T.^7 + ...
    1287.*(-1 + exp(T)).^5.*(1 + 247.*exp(T) + 4293.*exp(2.*T) + 15619.*exp(3.*T) + 15619.*exp(4.*T) + 4293.*exp(5.*T) + 247.*exp(6.*T) + exp(7.*T)).*t.^5.* ...
     T.^8 + 715.*(-1 + exp(T)).^4.*(1 + 502.*exp(T) + 14608.*exp(2.*T) + 88234.*exp(3.*T) + 156190.*exp(4.*T) + 88234.*exp(5.*T) + 14608.*exp(6.*T) + ...
      502.*exp(7.*T) + exp(8.*T)).*t.^4.*T.^9 + 286.*(-1 + exp(T)).^3.*(1 + 1013.*exp(T) + 47840.*exp(2.*T) + 455192.*exp(3.*T) + 1310354.*exp(4.*T) + ...
      1310354.*exp(5.*T) + 455192.*exp(6.*T) + 47840.*exp(7.*T) + 1013.*exp(8.*T) + exp(9.*T)).*t.^3.*T.^10 + ...
    78.*(-1 + exp(T)).^2.*(1 + 2036.*exp(T) + 152637.*exp(2.*T) + 2203488.*exp(3.*T) + 9738114.*exp(4.*T) + 15724248.*exp(5.*T) + 9738114.*exp(6.*T) + ...
      2203488.*exp(7.*T) + 152637.*exp(8.*T) + 2036.*exp(9.*T) + exp(10.*T)).*t.^2.*T.^11 + ...
    13.*(-1 + exp(T)).*(1 + 4083.*exp(T) + 478271.*exp(2.*T) + 10187685.*exp(3.*T) + 66318474.*exp(4.*T) + 162512286.*exp(5.*T) + 162512286.*exp(6.*T) + ...
      66318474.*exp(7.*T) + 10187685.*exp(8.*T) + 478271.*exp(9.*T) + 4083.*exp(10.*T) + exp(11.*T)).*t.*T.^12 + ...
    (1 + 8178.*exp(T) + 1479726.*exp(2.*T) + 45533450.*exp(3.*T) + 423281535.*exp(4.*T) + 1505621508.*exp(5.*T) + 2275172004.*exp(6.*T) + ...
    1505621508.*exp(7.*T) + 423281535.*exp(8.*T) + 45533450.*exp(9.*T) + 1479726.*exp(10.*T) + 8178.*exp(11.*T) + exp(12.*T)).*T.^13);
 g15=1./(-1 + exp(T)).^15.*exp(T).*((-1 + exp(T)).^14.*t.^14 + 14.*(-1 + exp(T)).^13.*t.^13.*T + 91.*(-1 + exp(T)).^12.*(1 + exp(T)).*t.^12.*T.^2 + ...
    364.*(-1 + exp(T)).^11.*(1 + 4.*exp(T) + exp(2.*T)).*t.^11.*T.^3 + 1001.*(-1 + exp(T)).^10.*(1 + 11.*exp(T) + 11.*exp(2.*T) + exp(3.*T)).*t.^10.*T.^4 + ...
    2002.*(-1 + exp(T)).^9.*(1 + 26.*exp(T) + 66.*exp(2.*T) + 26.*exp(3.*T) + exp(4.*T)).*t.^9.*T.^5 + ...
    3003.*(-1 + exp(T)).^8.*(1 + 57.*exp(T) + 302.*exp(2.*T) + 302.*exp(3.*T) + 57.*exp(4.*T) + exp(5.*T)).*t.^8.*T.^6 + ...
    3432.*(-1 + exp(T)).^7.*(1 + 120.*exp(T) + 1191.*exp(2.*T) + 2416.*exp(3.*T) + 1191.*exp(4.*T) + 120.*exp(5.*T) + exp(6.*T)).*t.^7.*T.^7 + ...
    3003.*(-1 + exp(T)).^6.*(1 + 247.*exp(T) + 4293.*exp(2.*T) + 15619.*exp(3.*T) + 15619.*exp(4.*T) + 4293.*exp(5.*T) + 247.*exp(6.*T) + exp(7.*T)).*t.^6.* ...
     T.^8 + 2002.*(-1 + exp(T)).^5.*(1 + 502.*exp(T) + 14608.*exp(2.*T) + 88234.*exp(3.*T) + 156190.*exp(4.*T) + 88234.*exp(5.*T) + 14608.*exp(6.*T) + ...
      502.*exp(7.*T) + exp(8.*T)).*t.^5.*T.^9 + 1001.*(-1 + exp(T)).^4.*(1 + 1013.*exp(T) + 47840.*exp(2.*T) + 455192.*exp(3.*T) + 1310354.*exp(4.*T) + ...
      1310354.*exp(5.*T) + 455192.*exp(6.*T) + 47840.*exp(7.*T) + 1013.*exp(8.*T) + exp(9.*T)).*t.^4.*T.^10 + ...
    364.*(-1 + exp(T)).^3.*(1 + 2036.*exp(T) + 152637.*exp(2.*T) + 2203488.*exp(3.*T) + 9738114.*exp(4.*T) + 15724248.*exp(5.*T) + 9738114.*exp(6.*T) + ...
      2203488.*exp(7.*T) + 152637.*exp(8.*T) + 2036.*exp(9.*T) + exp(10.*T)).*t.^3.*T.^11 + ...
    91.*(-1 + exp(T)).^2.*(1 + 4083.*exp(T) + 478271.*exp(2.*T) + 10187685.*exp(3.*T) + 66318474.*exp(4.*T) + 162512286.*exp(5.*T) + 162512286.*exp(6.*T) + ...
      66318474.*exp(7.*T) + 10187685.*exp(8.*T) + 478271.*exp(9.*T) + 4083.*exp(10.*T) + exp(11.*T)).*t.^2.*T.^12 + ...
    14.*(-1 + exp(T)).*(1 + 8178.*exp(T) + 1479726.*exp(2.*T) + 45533450.*exp(3.*T) + 423281535.*exp(4.*T) + 1505621508.*exp(5.*T) + ...
      2275172004.*exp(6.*T) + 1505621508.*exp(7.*T) + 423281535.*exp(8.*T) + 45533450.*exp(9.*T) + 1479726.*exp(10.*T) + 8178.*exp(11.*T) + ...
      exp(12.*T)).*t.*T.^13 + (1 + 16369.*exp(T) + 4537314.*exp(2.*T) + 198410786.*exp(3.*T) + 2571742175.*exp(4.*T) + 12843262863.*exp(5.*T) + ...
      27971176092.*exp(6.*T) + 27971176092.*exp(7.*T) + 12843262863.*exp(8.*T) + 2571742175.*exp(9.*T) + 198410786.*exp(10.*T) + ...
      4537314.*exp(11.*T) + 16369.*exp(12.*T) + exp(13.*T)).*T.^14);
   
   g1=repmat(g1,1,length(t));
   
   g=[g1' g2' g3' g4' g5' g6' g7' g8' g9' g10' g11' g12' g13' g14' g15'];
   