I have this error:
Error: File: p3.m Line: 118 Column: 1
Function definition not supported in this context.
Create functions in code file.
for this code:
clear all
clc

%% User Input

D = input('Diffusion coefficient (m^2/s): ');
t_end = input('End time (s): ');
L = input('Spatial dimension-L (m): ');
k = input('Reaction rate (1/s): ');
R = input('Retardation coefficient ®: ');
To = input('Tortuosity : ');
U = input('Velocity of channel flow : ');
boundary_condition = input('Boundary condition ([1]: No-Flux, [2]: Perfectly Absorbing): ');

%% Parameters

M = 1;
t_start = 1;
dt = 1;
dx = 0.1;
dy = 0.1;

time_range = t_start:dt:t_end;
x_range = [0:dx:2*L];
y_range = [-L:dy:L];

D_eff = D/R*To;
U_eff = U/R;

%% Initialize result matrix

Result = zeros(numel(y_range), numel(x_range), numel(time_range));

%% Compute concentrations

[X, Y, T] = meshgrid(x_range, y_range, time_range);
C_M = M ./ (4 * pi * T).^1.5 .* D_eff;


% Compute concentration based on the selected boundary condition

if boundary_condition == 1

    c = exp((-(X-U_eff*T).^2) ./ (4 * D_eff * T) + (-Y.^2) ./ (4 * D_eff * T) - k * T) + exp((-(X-U_eff*T).^2) ./ (4 * D_eff * T) + (-(Y + 2 * L).^2) ./ (4 * D_eff * T) - k * T) + exp((-(X-U_eff*T).^2) ./ (4 * D_eff * T) + (-(Y - 2 * L).^2) ./ (4 * D_eff * T) - k * T);

elseif boundary_condition == 2
    
    c = 0;
    for n = -1:1:1
    c = c + exp( (-(X - U_eff * T).^2)./(4 * D_eff * T) + (-(Y+ 4 * n *L).^2)./(4 * D_eff * T)) - exp( (-(X - U_eff * T).^2)./(4 * D_eff * T) + (-(Y+(4*n-2)*L).^2)./(4*D_eff*T));
    end

% Add continuous point source
c = c + f_source(X, Y, T);
    
%     c = 0;
%     for n = -1:1:1
%         c = c + exp( (-(X-U_eff*T).^2)./(4*D_eff*T) + (-(Y+4*n*L).^2)./(4*D_eff*T)) - exp( (-(X-U_eff*T).^2)./(4*D_eff*T) + (-(Y+(4*n-2)*L).^2)./(4*D_eff*T));
%     end

% else
% 
%     error('Invalid choice for boundary condition. Choose 1 for No-Flux or 2 for Perfectly Absorbing.');
% 
% end

Result = C_M .* c;

%% Visualization (2D Plot)
figure;
for t = time_range
    imagesc(Result(:, :, t));
    title(strcat('time: ', num2str(t), ' s'))

    xticks(1:(numel(x_range) - 1)/2:numel(x_range));
    xticklabels({num2str(-L), num2str(0), num2str(L)})
    xlabel('X-Axis (m)');

    yticks(1:(numel(y_range) - 1)/2:numel(y_range));
    yticklabels({num2str(L), num2str(0), num2str(-L)})
    ylabel('Y-Axis (m)');

    colorbar

    pause(0.01)
    colormap('jet');
end
close;

% 1D Concentration Profiles

x = round(numel(x_range) / 2);

t_1 = numel(time_range);
t_2 = round(numel(time_range) / 2);
t_3 = round(numel(time_range) / 4);

c_1 = Result(:, x, t_1);
c_2 = Result(:, x, t_2);
c_3 = Result(:, x, t_3);

figure;
plot(c_1, 1:numel(c_1));
hold on
plot(c_2, 1:numel(c_2)); 
hold on
plot(c_3, 1:numel(c_3));
hold off

xlim([0, 1.5 * max([max(c_1), max(c_2), max(c_3)])])
xlabel('Concentration (g/m^2)');
ylabel('Y-Axis (m)');
yticks(1:(numel(y_range) - 1)/4:numel(y_range));
yticklabels({num2str(-L), num2str(-L/2), num2str(0), num2str(L/2), num2str(L)})

legend({strcat('t = ', num2str(t_3), ' s'), strcat('t = ', num2str(t_2), ' s'), strcat('t = ', num2str(t_1), ' s') }, 'Location', 'northeast')

function c_source = f_source(X, Y, T)
    % Parameters for the continuous point source
    x_source = 0;                   % x-coordinate of the source
    y_source = 0;                   % y-coordinate of the source
    C_source = 1;                   % concentration of the source
    r_source = 0.5;                 % radius of the source

    % Calculate distance from each point to the source
    r = sqrt((X - x_source).^2 + (Y - y_source).^2);

    % Calculate concentration contribution from the source at each point
    c_source = (r <= r_source) .* C_source;   % Only points within the radius of the source will have a nonzero concentration
end