how to convert contour plot from mathematica to matlab

The thing to keep in mind is that Mathematica is formula-oriented and is programmed in terms of mathematical variables, while Matlab is numerically oriented and programmed in terms of matrices and arrays. To evaluate a function of two variables on a grid, it is necessary to first generate that grid explicitly in the form of arrays, and then to take care to use the element-wise operators (.*, ./ etc.) instead of the matrix operators. Here is the code:

[B, C] = meshgrid(b, c);

z = (-136*exp(power(88682,C) .* B) .* (1 - 1./exp(1*power(88682,C) .* B)))./(-1+ exp(power(88682,C) .* B)) ...
- 1*(power(3,C) + power(33,C) + power(146,C) + power(227,C) + power(342,C) + power(351,C) + power(353,C) + power(444,C) + power(556,C) + power(571,C) + power(709,C) + power(759,C) + power(836,C) + power(860,C) + power(968,C) + power(1056,C) + power(1726,C) + power(1846,C) + power(1872,C) + power(1986,C) + power(2311,C) + power(2366,C) + power(2608,C) + power(2676,C) + power(3098,C) + power(3278,C) + power(3288,C) + power(4434,C) + power(5034,C) + power(5049,C) + power(5085,C) + 2*power(5089,C) + power(5097,C) + power(2,1 + 3*C) .* power(5287,C) + power(5324,C) + power(5389,C) + power(5565,C) + power(5623,C) + power(6080,C) + power(6380,C) + power(6477,C) + power(6740,C) + power(7192,C) + power(7447,C) + power(7644,C) + power(7837,C) + power(7843,C) + power(7922,C) + power(8738,C) + power(10089,C) + power(10237,C) + power(10258,C) + power(10491,C) + power(10625,C) + power(10982,C) + power(11175,C) + power(11411,C) + power(11442,C) + power(11811,C) + 2*power(12559,C) + power(12791,C) + power(13121,C) + power(13486,C) + power(14708,C) + power(15251,C) + power(15261,C) + power(15277,C) + power(15806,C) + power(16185,C) + power(16229,C) + power(16358,C) + power(17168,C) + power(17458,C) + power(17758,C) + power(18287,C) + power(18568,C) + power(18728,C) + power(19556,C) + power(20567,C) + power(21012,C) + power(21308,C) + power(23063,C) + power(24127,C) + power(25910,C) + power(26770,C) + power(27753,C) + power(28460,C) + power(28493,C) + power(29361,C) + power(30085,C) + power(32408,C) + power(35338,C) + power(36799,C) + power(37642,C) + power(37654,C) + power(37915,C) + power(39715,C) + power(40580,C) + power(42015,C) + power(42045,C) + power(42188,C) + power(45406,C) + power(46653,C) + power(47596,C) + power(48296,C) + power(49171,C) + power(49416,C) + power(50145,C) + power(52042,C) + power(52489,C) + power(52875,C) + power(53321,C) + power(53443,C) + power(54433,C) + power(55381,C) + power(56463,C) + power(56485,C) + power(56560,C) + power(57042,C) + power(62551,C) + power(62651,C) + power(62661,C) + power(63732,C) + power(64103,C) + power(64893,C) + power(71043,C) + power(74364,C) + power(75409,C) + power(76057,C) + power(81542,C) + power(82702,C) + power(84566,C) + power(88682,C)) .* B ...
+ 1266.5546418385234*(-1 + C) + 136*log((136*exp(power(88682,C) .* B))./(-1 + exp(power(88682,C) .* B))) + 136*log(B) + 136*log(C);
contour(b, c, z)

I also took the liberty to slightly syntax-normalize and simplify your formula (e.g, power(2.718..., x) is exp(x)).

Here's the result:

If you want more contour lines, use something like contour(b, c, z, 100).

To understand what I've done here, read Array vs. Matrix operations and meshgrid from the Matlab documentation. Because .* means something else than * in Matlab, one should also take care not to have trailing decimal points in expressions, like 3.*B, though it doesn't really do any harm in that case.