MATLAB: Keep receiving error message when trying to create a surface plot of a a variable that is a function of two angles

errorsurface plot

I am analyzing the the forces observed on a rigid body that are a function of two angles, gamma and alpha. Using force and moment balances I created the following code which should solve for the desired variable (Gy) in terms of alpha and gamma. I then want to create a surface plot of Gy over a varying range of alpha and gamma angles but keep receiving error messages.
This is the code:
syms Gy; % substrate reaction force
syms Px Py; % shaft pin reaction forces
syms Wl Wc g; % weight of light and cartridge+clamp as well as gravity
syms S Sx Sy; % Spring forces
syms k x; % spring constant and spring displacement
syms gamma alpha; % angle with respect to ground (-45 to 45) and angle with respect to spring (0 to 25)
syms Dgpx Dgpy; % distance between point G and P for x and y axes
syms Dpsx Dpsy; % distance between point P and S for x and y axes
syms Dpwlx Dpwly; % distance between point P and Wl for x and y axes
syms Dpwcx Dpwcy; % distance between point P and Wc for x and y axes
syms b1 b2 b3; % variables for each value of the B vector in AX=B
Wl = 0.02; % In kilograms
Wc = 0.015; % In kilograms assuming material density of 1.25 g/cm^3 (PLA)
g = 9.81;
Dgpx = 15.3;
Dgpy = 12.7;
Dpsx = 0.8;
Dpsy = 27.8;
Dpwlx = 30;
Dpwly = 10.1;
Dpwcx = 13;
Dpwcy = 10.3;
x = 0.19846*(alpha).^(1.11259); % Using data points assembled curve of best fit for displacement as function of alpha
k = 0.4378; % spring model number 9434K23 - 2.5 lbs/in = 0.4378 N/mm
S = k*x; % Hooke's Law
Sx = S*cos(0.2617+alpha.*0.01745); % Spring force components

Sy = S*sin(0.2617+alpha.*0.01745); % Spring force components
% define matrix A: row 1 - sum of forces in y, row 2 - sum of forces in x,
% and row 3 - sum of moments about P
A = [1 1 0; 0 0 1; Dgpx 0 0];
% X = [Gy;Py;Px]
b1 = Sy+(Wl*g*cos(gamma*0.01745))+(Wc*g*cos(gamma*0.01745));
b2 = Sx-(Wl*g*sin(gamma*0.01745))-(Wc*g*sin(gamma*0.01745));
b3 = (Wl*g*(Dpwlx*cos(gamma*0.01745)-Dpwly*sin(gamma*0.01745)))+(Wc*g*(Dpwcx*cos(gamma*0.01745)-Dpwcy*sin(gamma*0.01745)))+(Sx*Dpsy)-(Sy*Dpsx);
B = [b1;b2;b3]; % assign B vector values
X = linsolve(A,B); % solve AX=B for X
Gy = X(1,1);
Py = X(2,1);
Px = X(3,1);
SubForce(1:19,1:6) = Gy;
[gamma,alpha] = meshgrid(-45:5:45,0:5:25);
surf(gamma,alpha,SubForce);
Any help would be greatly appreciated. I am not good at using Matlab at all and can't seem to figure out what the issue is.
Thanks.

Best Answer

Change
SubForce(1:19,1:6) = Gy;
[gamma,alpha] = meshgrid(-45:5:45,0:5:25);
surf(gamma,alpha,SubForce);
to
[Gamma, Alpha] = meshgrid(-45:5:45,0:5:25);
SubForce = double(subs(Gy, {gamma, alpha}, {Gamma*pi/180, Alpha})); %angles need to be radians for cos and sin
surf(Gamma, Alpha, SubForce)