How to substitute the value in symbolic function?

I have a problem with MATLAB coding. I'm trying to find the solution for the fractional-order equation. I already found for c0 until c3. Then I find sum of series as follow, and m=0,1,...,4. The problem is I dont know how to substitute the value in symbolic function and the command window shows 'Inputs must be a scalar and a square matrix' . I dont really understand what it's mean. Why this happened? And how to fix my code? I hope someone can help me.

% c(x,t) is the concentration of GBM tumour cells
clear;

syms x t D n h q

% Parameters to define the diffusion equation and the range in space and
% time
a = 1;  % Step size for time
b = 0.5;    % Step size for place
t = 0:a:1500;   % Time points
N = length(t);  % Length of time points
r = 0:b:50; % Place points
R = length(r);  % Length of place points

% Parameters needed to solve the equation
p = 0.012;  % Proliferation rate (per day)
x0 = 25;    % The middle of considered interval (mm)
epsilon = 0.01; % suggest by Ozugurlu
cmax = 62.5;    % Carrying capacity of GBM tumour cells
alpha = 1.0;    % Fractional derivative order
q = 1/n;

M = input('Enter length of series (M) =');
sum_c1 = 0;

% Initial concentration
    c(1) = (1/((sqrt(2*pi()))*epsilon))*exp((-1/2)*(((x-x0)/epsilon)^2));   % c0
    c(2) = (D/((sqrt(2*pi()))*(epsilon^3))-(D*((x0-x)^2))/((sqrt(2*pi()))*(epsilon^5))-p/((sqrt(2*pi()))*epsilon))*(h*((t^alpha)/(gamma(alpha+1))))*(exp((-1/2)*(((x-x0)/epsilon)^2)))+((h*p*(t^alpha))/(cmax*2*pi()*(epsilon^2)*(gamma(alpha+1))))*exp(-1*(((x-x0)/epsilon)^2));    % c1

for m=2:M
    
    derivative = diff(c(m),x,2);
    c(m+1) = (n+h)*c(m)-(h*(t^alpha)/(gamma(alpha+1)))*(D*derivative+p*c(m)-(p/cmax)*((c(m))^2))
    c1(m+1) = c(m+1)*(q^m)
    
end

sum_c1 = sum_c1+c1(m+1)
cn = c(1)+c(2)*q+sum_c1

for j=1:N
    
    t(j+1) = t(j)+a;
    
    for i=1:R
        
        x(i+1) = x(i)+b;
        
        if i<=16
            D(i+1) = 0.13;
        elseif 16<=i && i<=86
            D(i+1) = 0.65;
        else
            D(i+1) = 0.13;
        end
    end
    
    cn(i+1,j+1) = cn;
    
    subs(cn,{x,t,D,n,h},{x(i+1),t(j+1),D(i+1),1,-1});
    double(subs(cn,{x,t,D,n,h,q},{x(i+1),t(j+1),D(i+1),1,-1}));
end

[KSSV]

Read about subs.