As already pointed out by Walter Roberson, 4 points cannot be coplanar. Assuming that they are, you may try this
clear variables, close all
P = rand(3,3);
V1 = P(2,:)-P(1,:);
V2 = P(3,:)-P(1,:);
P = [P; P(1,:)+V1+V2];
figure, hold on, axis equal
plot3(P(:,1),P(:,2),P(:,3),'o');
quiver3(P(1,1),P(1,2),P(1,3),V1(1),V1(2),V1(3),'AutoScale','off')
quiver3(P(1,1),P(1,2),P(1,3),V2(1),V2(2),V2(3),'AutoScale','off')
B0 = mean(P,1);
plot3(B0(1),B0(2),B0(3),'*');
W = cross(V2,V1);
U = W/sqrt(sum(W.^2,2));
k = 2;
B = B0+k*U;
quiver3(B0(1),B0(2),B0(3),k*U(1),k*U(2),k*U(3),'AutoScale','off');
plot3(B(1),B(2),B(3),'s');
Basically two vectors are created from the first three points, assuming the first point as reference. The orthogonal direction is obtained with a cross product of these vectors. Note that you can control the positive side by changing the order of the cross product.
Best Answer