how to calculate the drivative of discretized ODE

I am discretizing DDE to ODE using pseudospectral method. I want to compute derivative of its solution for training state and want to use the right hand equation of the discretized ODE here dODE represents discretized ODE

tspan=[0 20];  M=5; t_r=0.8;
phi = @(x) cos(x);  
g = @(t,y,Z,par) par * y * (1 - Z);
tau = 1;
par = 1.5;
[D, theta] = difmat(-tau, 0, M);
X = phi(theta); 
u0=X;
options = odeset('RelTol', 1e-10, 'AbsTol', 1e-10);
sol = ode45(@(t,u) dODE(t, u, par, tau, M,D), tspan, u0,options);
x_n = sol.y ;
t = linspace(tspan(1), tspan(2), 100);
x_an = interp1(sol.x, x_n(1,:), t, 'linear');
n_s_r = size(x_an, 2);  
x_a = floor(n_s_r * t_r);
x_tn = x_an(:, 1:x_a);
n_all = length(t);
n_train = round(t_r * n_all);
t_t = t(1:n_train);
function dydt = dODE(t,u, par, tau, M, D)
yM = u(1);
VM = u(2:end);
dMDM_DDE = kron(D(2:end,:), eye(1));
dydt = [par*yM*(1-VM(end)); (dMDM_DDE)*[yM;VM]];
end
function [D, x] = difmat(a, b, M)
% CHEB pseudospectral differentiation matrix on Chebyshev nodes.
% [D,x]=CHEB(a,b,M) returns the pseudospectral differentiation matrix D

if M == 0
    x = 1;
    D = 0;
    return
end
x = ((b - a) * cos(pi * (0:M)' / M) + b + a) / 2;
c = [2; ones(M-1, 1); 2].*(-1).^(0:M)';
X = repmat(x, 1, M+1);
dX = X - X';
D = (c * (1./c)')./(dX + (eye(M+1)));
D = D - diag(sum(D'));
end


(like this way 
for j = 1:length(t_t)
        DX(:,j) = g(t_t(j), x_tn(j), x_dn(:,j), par);
end this is derivative of original DDE)

Answer 1

[Torsten]

After the line

sol = ode45(@(t,u) dODE(t, u, par, tau, M,D), tspan, u0,options);

you can compute the derivatives of the solution as

[~,yp] = deval(sol,sol.x) 

[Muhammad]

Thank you for reply, but i want to comupte the derivative of this right hand side

dODE(t,u, par, tau, M, D).

[Torsten]

Here is an example:

fun = @(t,y) [y(1);-y(2)];
tspan = 0:0.1:1;
y0 = [1;1];
sol = ode45(fun,tspan,y0);
% Compute 1st derivative of the solution
[~,yp] = deval(sol,sol.x);
% Compute 2nd derivatve of the solution
n = size(yp,2);
ypp(:,1)=(yp(:,2) - yp(:,1))./(sol.x(2)-sol.x(1));
ypp(:,2:n-1) = (yp(:,3:n) - yp(:,1:n-2))./(sol.x(3:n)-sol.x(1:n-2));
ypp(:,n)=(yp(:,n) - yp(:,n-1))/(sol.x(n) - sol.x(n-1));
% Plot the 2nd derivative of the solution
plot(sol.x,ypp)
grid on

[Muhammad]

Thank you for providing this