Here is MATLAB's meshgrid. I came up with these implementations in Mathematica:

{x, y} = With[{n = 10}, Table[#, {i, n}, {j, n}]] & /@ {j, i};


{x, y, z} = Array[List, {5, 5, 5}][[All, All, #]] & /@ {3, 2, 1}

but it seems that it's not fast enough. The following example is about three times slower than MATLAB.

(* Mathematica *)
  {a, b} = With[{n = 3000}, Table[#, {i, n}, {j, n}]] & /@ {j, i};
  c = a^2 + b^2 // N // Sqrt;
  Compile[{}, Position[Unitize@FractionalPart@c, 0]][] // Length
) // AbsoluteTiming

Can you recommend an efficient method?


  {a, b} = {#, Transpose@#} &@ConstantArray[Range@3000, {3000}];
  c = Sqrt@N[a*a + b*b];
  SparseArray[Unitize@FractionalPart@c~BitXor~1]["NonzeroPositions"] // Length
) // AbsoluteTiming
  • 8
    $\begingroup$ For those of us who do not use Matlab can you describe what the meshgrid function does ...rather than forcing us to infer it by examining your code. $\endgroup$ Commented Aug 19, 2013 at 4:41

4 Answers 4


(Update, added more points, and more timings)

Using MATLAB's help standard example for meshgrid:

Mathematica implementation

meshgrid[x_List, y_List]:={ConstantArray[x,Length[x]],Transpose@ConstantArray[y,Length[y]]}

{xx, yy} = meshgrid[Range[-2, 2, .1], Range[-4, 4, .2]];
c = xx*Exp[-xx^2 - yy^2];

pts = Flatten[{xx, yy, c}, {2, 3}];
ListPlot3D[pts, PlotRange -> All, AxesLabel -> Automatic, 
           ImagePadding -> 20, Mesh -> 35, InterpolationOrder -> 2, 
           ColorFunction -> "Rainbow", Boxed -> False]


    [X,Y] = meshgrid(-2:.1:2, -4:.2:4);
    Z = X .* exp(-X.^2 - Y.^2);

Enter image description here

MATLAB timing was done using this template code by changing the grid spacing, as an example for the code that generate the above plot:

tic;  [X,Y] = meshgrid(-2:.2:2, -4:.4:4); toc

Mathematica was done on this template

 {xx, yy} = meshgrid[Range[-2, 2, .2], Range[-4, 4, .4]];

Table of timings

The timing were done for different grid sizes as in (.2,.4) pairs in the above example and using Mathematica 9.01 and MATLAB 2012a, all on Windows 7. I could not do smaller grids than shown in the table below, else I ran out of memory. All times are in seconds. Same PC. 16 GB RAM, Intel Core i7

Mathematica graphics

Mathematica's above implementation of meshgrid using ConstantArray seems to be a little slower as the number of grid points increased. A faster implementation should be investigated.

Mathematica graphics

  • 1
    $\begingroup$ I think ListPlot3D[Transpose[c],...] is enough. $\endgroup$
    – user21
    Commented Aug 19, 2013 at 8:52
  • 1
    $\begingroup$ You could also use pts = Flatten[{xx, yy, c}, {2, 3}] $\endgroup$ Commented Aug 20, 2013 at 10:34
  • 3
    $\begingroup$ @SimonWoods thanks! that is even simpler. Just updated. I must say I find Matlab little easier when it comes to this part. Matlab wants each coordinates in one separate vector, so one passes 3 separate vectors, while Mathematica wants the coordinates in a matrix (each row has each {x,y,z} coordinate). I never understood why can't there be 2 API's for this in Mathematica. Where the user can choose to pass the coordinates as {X...},{Y...},{Z..} vs {{X,Y,Z}...}, and let the user choose the simple one depending. What is so hard about having two different interfaces like this for the command? $\endgroup$
    – Nasser
    Commented Aug 20, 2013 at 14:56
  • 1
    $\begingroup$ You can also use ListPlot3D[Transpose@c, DataRange -> {{-4, 4}, {-2, 2}}, PlotRange -> All] $\endgroup$
    – chyanog
    Commented Sep 2, 2013 at 2:29
  • 1
    $\begingroup$ At large n I think it's more efficient to build the list via Map than to Transpose it. i.e. use Map[ConstantArray[#, Length[y]] &, y] instead of Transpose@ConstantArray[y, Length@y]. It's not much better, but it does give a second or so of improvement when you have 15000 mesh points in both x and y $\endgroup$
    – b3m2a1
    Commented Feb 3, 2018 at 0:11


Another way to define meshgrid() in Mathematica is:

meshgrid[xgrid_List, ygrid_List] := Transpose[Outer[List, xgrid, ygrid], {3, 2, 1}]

It's definitely not winning the speed competition (about 10x slower than rm -rf's solution), but speaks more for Mathematica's language, which allows neat and concise definitions like this.


Regarding speed, in your particular example you can do the squaring before the repetition via ConstantArray, loose the superfluous N in the computation of c and use Count instead of Position[...]\\Length to save computation time.

({a, b} = {#, Transpose@#} &@ConstantArray[N@Range@3000, {3000}];
 c = a^2 + b^2 // N // Sqrt;
 Compile[{}, Position[Unitize@FractionalPart@c, 0]][] //Length) // AbsoluteTiming
(* {0.698965,7388} *)

({asq, bsq} = {#, Transpose@#} &@ ConstantArray[N[Range[3000]]^2, {3000}];
 c = Sqrt[asq + bsq];
 Compile[{}, Count[Unitize@FractionalPart@c, 0, 2]][]) // AbsoluteTiming
(* {0.396491,7388} *)
  • $\begingroup$ (+1) I contemplated posting the same thing using Outer (which IMO is the clearest way to think about tensored grids) but I noticed that this turned into a speed chase and got discouraged. $\endgroup$
    – gpap
    Commented Aug 20, 2013 at 16:42

The accepted answer excellent and comprehensive, but it has a very fine/little mistake\bug. I've tried to fix it, but it has to be accepted due to my current reputation(=101) on this site(well, I'm not very active here), and mistakenly it was not. The mistake is in the correct treatment of lists of different sizes. I'll explain in more details...

consider two implementations, the first (meshgrid1) is copied from current version (as of 2/2/2018) of the accepted answer, the second (meshgrid2) is a correctly reproduces matlab's meshgrid. The little difference is in the argument of Length.

meshgrid1[x_List, y_List]:={ConstantArray[x,Length[x]],Transpose@ConstantArray[y,Length[y]]}

meshgrid2[x_List, y_List]:={ConstantArray[x,Length[y]],Transpose@ConstantArray[y,Length[x]]}

now use it in a slightly different example then in the accepted answer (the ylist here is not the same size as xlist)

xlist = Range[-2, 2, .1];   
ylist = Range[-4, 4, .1];   

Out: 41
Out: 81

and see the dimensions of the result

{xx1, yy1} = meshgrid1[xlist, ylist];   

Out: {41, 41} 
Out: {81, 81}


{xx2, yy2} = meshgrid2[xlist, ylist];

Out: {81, 41}
Out: {81, 41}

now the same thing in matlab - decide by yourself which one is correct or make more sense

>> xlist = -2:.1:2; 
>> ylist = -4:.1:4;
>> [xx,yy] = meshgrid(xlist, ylist);

>> size(xlist)

ans =

     1    41

>> size(ylist)

ans =

     1    81

ans =

    81    41

>> size(yy)

ans =

    81    41
  • 1
    $\begingroup$ Good job on debugging. The accepted answer probably got confused since the two starting vectors happened to be the same length, so the wrong lengths being supplied was not noticed. Had a rectangular plotting region been used like in your case, the error would have been spotted more quickly. $\endgroup$ Commented Mar 21, 2018 at 4:59
  • 1
    $\begingroup$ This answer should be accepted and be seen by more people who really want to use the meshgrid fucntions. $\endgroup$
    – Knife Lee
    Commented Feb 12, 2022 at 11:11

meshgrid()'s functionality is easily handled by ConstantArray. For the example you provided (i.e. meshgrid called with a single vector),

{a, b} = {#, Transpose@#}&@ConstantArray[N@Range@3000, {3000}]; // AbsoluteTiming
(* {0.068030, Null} *)

You can easily extend it to handle the two argument case by calling ConstantArray for each list (x and y) and using the length of the other list as the number of repetitions, which I'll leave for you to implement.

  • 1
    $\begingroup$ This brings down the time from ~0.9s to ~0.5s on my system. Still slower than MATLAB... I think an explicit Position finder (i.e., using Do loops) that is compiled to C in Mathematica might speed things up a bit more... $\endgroup$
    – rm -rf
    Commented Aug 19, 2013 at 4:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.