Here is code to find the unit normal to a plane defined by three points, and then the angle with the z-axis.. In the example below, the points are chosen at random. The gamma=min(..) statement assures that the answer will be in the range 0 <= gamma <= pi/2.
p0=rand(3,1); p1=rand(3,1); p2=rand(3,1);
fprintf('Angle with Z-axis=%.3f radians\n',gamma)
I asked whether you know the x,y,z coordinates exactly, because the answer affects you choice of what to minimize to get the "best fit". A. If you know the x and y values exactly, andt there is measurement error or of other uncertainty along z, then the best-fit is the fit that minimizes the sum squared error in the z-direction. This is ordinary least squares, or standard multiple linear regression. B. If you have measurement errors or other uncertinty on all three axes, then the best fit is the fit that minimizes the orthogonal sum squared error from the plane, also known as orthogonal or total least squares fitting.
A. Ordinary least squares=standard multiple regresion: If you have access to the Statistics toolbox, use fitlm(). If you don;t, you can do it the old fashioned way. I will put code for both in a separate post.
B. Orthogonal (total) least squares: use pca() to find the principal components (PCs). There will be thee PCs. The lthird (last) PC, which accounts for the smallest portion of the variance, is the normal to the best fit plane. I will put code in a separate post.
Once you have a unit normal vector n, from method A or B, then the angle with the Z axis is given as before:
Code for A and B in separate post.