Greetings --
Apologies if this is posted to the wrong forum. Have been using Maple and Maxima for 20+ years, and have only recently started using Mathematica (v. 14). Having a heck of a time getting Mathematica to do what is trivial in either Maple or Maxima (not a slam on Mathematica - just an honest statement of where I am on the learning curve). Simple example: 2 equations in 2 unknowns (say, x1 and x2). Want to derive the Jacobian of these equations wrt x1 and x2, and print said Jacobian in matrix form. This is basically 3-4 lines in Maple or Maxima. But in Mathematica? Here is what I've tried.
(*Define the functions*)
f1[x1_, x2_] := x1^2 + x2^2
f2[x1_, x2_] := x1 x2
(*Compute the Jacobian matrix with respect to x1 and x2*)
jacobian = JacobianMatrix[{f1[x1, x2], f2[x1, x2]}, {x1, x2}]
(*Display the Jacobian matrix*)
MatrixForm[jacobian]
But, the final command to display the Jacobian returns nada of any use. It simply returns the following:
JacobianMatrix[{x1^2+x2^2,x1x2},{x1,x2}]
So, how the heck do I get Mathematica to output what Maple (for example - same code more or less in Maxima) outputs in one command, looking like what I expect (i.e., the formatted Jacobian matrix, below):

Many thanks in advance. I was warned the learning curve for Mathematica was appreciably steeper than (say) Maple, but I wasn't expected even something this simple to stump me quite so soon.