Derivative after parameter - system of ode \begin{align}
\frac{dx_{1}}{dt}& = x_{2}
\\
\frac{dx_{2}}{dt} &= - A x_{1} - B x_{2} + a\sin(Ct)
\end{align}
$x_{1}(0)=1$ and $x_{2}(0)=0$ 
I want to draw 2 plots in Matlab: $\frac{dx_{1}}{dA}$ and $\frac{dx_{2}}{dA}$
I found: $\frac{dx_{1}}{dA}$ = t 
and $\frac{dx_{2}}{dA}$ = $\int_{0}^{t} x_{2} e^t ds$ 
How can i solve $\frac{dx_{2}}{dA}$ ?
Can somebody tell me how i can do it in matlab?
 A: Whenever you have a parameter depending ODE
$$
\dot x=f(t,x,q),~~~x(0,q)=x_0(q)
$$
then the derivation of the solution for the parameters takes the form
$$
\frac∂{∂t}\frac∂{∂q_k}x(t,q)=\frac∂{∂x}f(t,x,q)\frac∂{∂q_k}x(t,q)+\frac∂{∂q_k}f(t,x,q),
\\~\\
\frac∂{∂q_k}x(0,q)=\frac∂{∂q_k}x_0(q),
$$
or perhaps more readable, set $u_k=\frac∂{∂q_k}x$ and use "operator" notation, then
$$
∂_tu_k(t,q)=∂_xf(t,x,q)u_k(t,q)+∂_{q_k}f(t,x,q),
\\~\\
u_k(0,q)=∂_{q_k}x_0(q).
$$
Note that in general this equation contains the original solution $x$, so it increases the system. In a numerical solver, the new equations can be solved along the original equations with a correspondingly increased state space vector.

In the current case, set $u=\frac{∂x}{∂A}$, then
\begin{align}
\dot u_1&=u_2,&~~u_1(0)=0,\\
\dot u_2&=-x_1-Au_1-Bu_2& u_2(0)&=0.
\end{align}
This in no way has the stated solutions. In matlab you would formulate the extended system like
function xdot=ext_sys(t,x,a,A,B,C)
  u = x(3:4);
  dx = [ x(2); -A*x(1)-B*x(2)+a*sin(C*t) ]
  du = [ u(2); -x(1)-A*u(1)-B*u(2) ]
  xdot = [ dx ; du ]


(Update 2/27/2021) Debugging Kutta's 3/8 method in the comments with better (slightly) formatting
function [t, x] = RK38th(f, x0, t0, T, h)
  t = t0:h:T;
  nt = numel(t);
  nx = numel(x0);
  x = nan(nx, nt);
  x(:,1) = x0;
  for k = 1:nt-1
    k1 = h*f(t(k), x(:,k));
    k2 = h*f(t(k) + h/3 ,  x(:,k) + k1/3 );
    k3 = h*f(t(k) + (2*h)/3, x(:,k) - k1/3 + k2);
    k4 = h*f(t(k) + h, x(:,k) + k1 - k2 + k3);
    dx = (k1 + 3*k2 + 3*k3 + k4)/8;
    x(:,k+1) = x(:,k) + dx;
  end
end

Note that in k2, the x update is now like in the other stages composed of the k vectors.
This then results in the following plot, where the left part shows the solution, and the right part compares the computed derivative with a difference quotient, where the difference in $A$ was chosen large as $0.1$ to get visible deviations. For smaller values of eps, the curves are visually the same.

The parameters and calling code are
  A = 2; B = 0.02; C = 1; a = 0.8;
  t0 = 0; T = 10; h = 0.5;

  x0 = [2; 0];
  x0_ext = [ x0(1); x0(2); 0; 0 ];
  
  eps = 1e-1;
  [T0,X0] = RK38th(@(t,x)ext_sys(t,x,a,A,B,C), x0_ext, t0, T, h);
  [T1,X1] = RK38th(@(t,x)ext_sys(t,x,a,A+eps,B,C), x0_ext, t0, T, h);

  subplot(121);
  plot(T0,X0);
  legend(['x1'; 'x2'; 'u1'; 'u2']);

  subplot(122);
  plot(T0,(X1(1:2,:)-X0(1:2,:))/eps, '-+', T0, X0(3:4,:));
  legend(['dx1'; 'dx2'; 'u1'; 'u2']);

