Please check me.May I show you two figures.The first figure is from my matlab code. But I want to obtain second figure at time t=t1, x-axis from -50 to 150, y-axis form -150 to 150 and z-axis form 0 to 0.5 for one wave. But I didn't get.

%%% Illustration of the free surface elevation %%%
clear all
close all
clc
zeta0 = 1;     % Maximum displacement(m)
L = 100;      % Propagate length (km)
W = 100;      % Source width (km)
v = 8.4;     % Ruptuer velocity to obtain maximum surface amplitude(km/m)
h = 2;        % Water depth (km)
g = 0.0098;   % Acceleration due to gravity (km/s^2)
% g = 35.28;   % Acceleration due to gravity (km/m^2)
% t1 = 50/v;    % (min)
% t = t1;     % (min)

%% Free surface elevation %%
syms X Y k1 k2 s t
T = ((zeta0*v*t)/(2*L))*(1-cos((pi.*X)/50)).*(1-(cos(pi.*(Y+150))/100));
T1 = zeta0*v*t/(2*L)*(1-cos(pi/50*X)).*(1-cos(pi/100*(Y+150))) ; 
T2 = zeta0*v*t/L*(1-cos(pi/50*X)) ;
T3 = zeta0*v*t/(2*L)*(1-cos(pi/50*X)).*(1-cos(pi/100*(Y-150))) ; 
T = T1+T2+T3
% Laplace transform of 'T'.
zeta_1 = laplace(T,t,s);

% Double Fourier transform of 'eta_1' for var 'X' and 'Y' to transVar 'k1' and 'k2'.
zeta_2 = fourier(zeta_1,X,k1);
zeta_ = fourier(zeta_2,Y,k2);

k = sqrt(k1^2+k2^2);
w = sqrt(g*k*tanh(k*h));

% Substituting 'eta_' into 'zeta_1' formula.
eta_1 = (zeta_.*s^2)/(cosh(k*h)*(s^2+w^2));

% Inverse Laplace transform of 'zeta_1' for var 's' to transVar 't1'.
eta_ = ilaplace(eta_1,s,t);

% Double inverse Fourier transform of 'zeta_' for var 'k1' and 'k2' to transVar 'X' and 'Y'.
eta1=ifourier(eta_,k1,X);
eta=ifourier(eta1,k2,Y);

%% For figure %%
Xvals = -50:150;
Yvals = linspace(-150,150,5); 
[X1,Y1] = meshgrid(Xvals,Yvals);
Zsym = subs(eta, {X, Y, t}, {X1, Y1, 50/v});
Z = double(Zsym);
surf(X1, Y1, abs(Z))
% axis([-200 200 -200 200 0 3])
xlabel('x(km)','FontSize',15)
ylabel('y(km)','FontSize',15)
title('t = t1 min','FontSize',15)