Linear regression with multiple variables and multiple unknown coefficients

I have this formula, where is the Boltzmann's constant, is the melting point, b and W are the base and width of a 2D interface, respectively. I have a data file that produces x, y, and z values that corrospond to , , and , respectively as well. I need to determine the value of the coefficients and . I tried to do a linear regression for and separately, but I don't believe this is the right approach. Is there a way to do a multiple varable linear regression to solve for multiple coefficients in MATLAB?

To describe the code a bit, a loglog plot was made of and . The solid line represents data points that follow . The data values between the vertical dashed lines are the only ones needed to be studied; everything else can be ignored.

fourier = readtable('(110).fourier.dat');
fourier = fourier{:,:};

A_k = reshape(fourier(:,3),62,62);
A_kx = mean(A_k); A_ky = mean(A_k,2);
k_b = 1.38064852e-23;                               %Boltzmann's constant: J/K
T_m = 957;                                          %Melting point temperature: K
b = 60.385;                                         %Base: Å; x-axis cell size
W = 61.1103;                                        %Length: Å; y-axis cell size
k_x = fourier(1:62,1)';k_y = fourier(1:62:end,2);

subplot(1,2,1)
plot(log10(k_x(2:end)),log10(A_kx(2:end)),'o','linewidth',.75), hold on
xlabel('log_{10}k_x    [Å^{-1}]'), ylabel('log_{10}< |A(k)|^2 >    [Å^2]')

s = 5; t = 15;
xline(log10(k_x(s)),'k--','linewidth',3),xline(log10(k_x(t)),'k--','linewidth',3)

The next bit of code was trying to find a linear regression of the selected data points. The slope of this regression line was believed to equal .

y = k_b*T_m/(L_x*L_y)*A_kx.^(-1);

subplot(1,2,2)
p = k_x(s:t).^2;
q = y(s:t);
g = polyfit(p,q,1);
ghat = [g(1),g(2)];
h = polyval(ghat,p);

plot(p,q,'ok','linewidth',0.5,'MarkerFaceColor',[0 0 0]), hold on
plot(p,h,'-k','linewidth',1.5), hold off
axis square, set(gca,'FontSize',18,'LineWidth',1,'TickLength',[0.025 0.025])
xlabel('k_x^2    [Å^{-2}]'); ylabel('k_bT_m / bW <|A(k)|^2>    [J/Å^2]')

I can do this for both the x and the y, but again, I don't think this is correct. Is there a way I can do both at the same time in MATLAB?