How to Replace Derivative Terms in Equations

clc,clear
%已知氢气在空气中浓度达到4.0%～75.6%就会爆炸
%假设初始条件：t = 0; H2 : O2 : N2 = 0.1 : 0.198 : 0.702
c = sym('c',[3,1]);%the units of concentration should be volume fraction
syms t z
c(1) = str2sym('c1(t,z)');%h2
c(2) = str2sym('c2(t,z)');%o2
c(3) = str2sym('c3(t,z)');%n2
c0 = [0.1;0.198;0.702];
% c = [c1;c2;c3];
%%
%parameters
%units of mass:kg; units of distance:km;
%moleculer weight
m_H = 1.67E-27;%kg
m = [2 * m_H;32 * m_H;28 * m_H];
%reduced mass
for i = 1:3
    for j = 1:3
        if i~=j
            mu(i,j) = m(i) * m(j) / (m(i) + m(j));
        else
            mu(i,j) = 0;
        end
    end
end
%%
%the sum of kinetic raddi
k_r = [289E-15;364E-15;346E-15];%km
for i = 1:3
    for j = 1:3
        if i~=j
            sigma(i,j) = k_r(i) + k_r(j);
        else
            sigma(i,j) = 0;
        end
    end
end
%%
%Boltzmann constant
k = 1.381E-23;%J/K
%T-z
dTdz = -6.5;%K/km
T0 = 293.15;%K
T = T0 + dTdz * z;
%scale height
g = 9.8;%m/s^2
R = 8.314;%J/mol/K
M = [2;32;28];
sum_c = sum(c);
sum_c0 = sum(c0);
%%
for i = 1:3
    %thermal diffusion factor
    alpha(i) = -log((c(i) ./ c0(i)) .* ((sum_c0 - c0(i)) ./ (sum_c - c(i)))) ./ log(T ./ T0);
    %individual height
    Hi(i) = (R * T) / (M(i) * g);
    for j = 1:3
        if j ~= i
            D_d(i,j) = c(j) .* (16/3) .* (mu(i,j)./m(i)) .* sigma(i,j).^2 .* (pi.*k.*T./(2.*mu(i,j))).^0.5
        end
    end
end
for i = 1:3
    % moleculer diffusion coefficent
    sumD_d(i) = sum(D_d(i,:))
    D(i) = Hi(i) .* g ./ sumD_d(i)
    %velocity-moleculer diffusion
    v_m(i) = D(i) .* ((1 ./ c(i))*diff(c(i),z) + (1 + alpha(i)) ./ T .* dTdz + 1 ./ Hi(i))
end
for i = 1:3
    %eddy diffsion coefficient--< Transport phenomena and the chemical composition  in the mesosphere and lower thermosphere>
    K_M = 3.8E6;%cm2s-1
    K_m = 6E5;
    S2 = 0.18;%km-2
    S3 = 0.098;%km-1
    z0 = 102;%km
    K = (K_M - K_m) * exp(-S2 * (z - z0).^2) + K_m * exp(-S3 * (z - z0));
    M_air = 29;
    H = (R * T) / (M_air * g);
    %velocity-eddy diffusion
    v_e(i) = K .* ((1 ./ c(i)) .* diff(c(i),z) + dTdz ./ T + 1 ./ H)
end
v = -v_m - v_e;
%%
for i = 1:3
    eqn(i) = diff(c(i),t) == -diff(c(i) .* v(i),z)
end
%%
syms delta_z
C = [];
% for i = 1:3
    eqn1 = subs(eqn(1),diff(c1(t, z), z),(C(i+1,n) - C(i,n)) ./ delta_z)
    

Matlab discretes the differential equation and replaces the differential term with the differential term. How can we do it ?The subs is not useful.

That is in an equation, replace the differential term on the left with the difference term on the right.

[KSSV]

Read about diff.

[Torsten]

What do you intend with your symbolic setup ?

If in the end, you replace the differentials with finite-difference approximations and calculate the concentrations numerically, I'd skip all the symbolic calculations and set up the problem numerically right from the beginning.

[Hou X.Y]

Yeah,I want to use finite-difference approximations to calculate the concentrations numerically,and the initial equation is complicated,so I use the symbolic calculations to deduce the final term of the diffrential function which can be used for calclations numerically.