I made some plots for 1d Ising chain with finite N, and it seems like there is always a maximum of specific heat and susceptibility at certain temperature. As the N gets larger and larger, the maximum moves to T=0. We know there is no phase transition in 1d Ising model, then what the temperatures of maximum of specific heat and susceptibility correspond to? I attach some plots below for 4 spins Ising chain with periodic boundary condition.
This behavior is expected, and it's due to the finite-size scaling. You will see the same behavior for 2D case as well. As you make the system larger, your results will get closer to what you expect from analytical solution. After all, analytical solutions are derived for infinite chains/lattices. No wonder if a small scale lattice can't simulate their behavior completely. And because of this, we tend to use Monte Carlo or other methods for large enough lattices.
check this for example: http://physics.bu.edu/~py502/slides/l17.pdf