IRFZ44N — хороший транзистор, но может быть ему недостаточно вашего напряжения, а 3.3 вольта уж точно мало. Возьмите какой-нибудь из серии IRL, я раньше часто использовал IRL2703 или IRL2705.
Не забывайте смотреть на график №3 в даташитах, да и вообще на сам даташит.

У вас сейчас много точек отказа, их нужно по отдельности проверить. Попробуйте включить/выключить транзистор сами, подавая на затвор 0 или 10 вольт. Затвор при этом должен быть замкнут на землю через резистор больше 10кОм, чтобы накопленный заряд рассасывался.
Проверить транзистор нужно обязательно, потому что они довольно чувствительны к статике, и обижаются когда их трогают мокрыми руками.

Я не использую ардуино (более того, я ардуинохейтер), но "тривиальная смена analogWrite(pinNumber, HIGH) и analogWrite(pinNumber, LOW) с частотой в 1000 мс" — разве это правильно? По-моему, в этом варианте нужно использовать digitalWrite.

Дурацкая какая-то схема, что значит "0-60V+" и "0-60V-"? Это явно от недомыслия автора схемы. Вы не могли из-за этого случайно перепутать + и -?

P.S. почитал сайт ergoz.ru, это просто насилие над мозгом. Нельзя обучаться по таким руководствам, разве что если себя не уважать, весь текст — дословный перевод с английского:
* изобретение самодельных терминов, например "стягивающий резистор" или "ворота";
* сплошные кальки с английского (например, "Arduino действительно пульсирует от 0 до 5V");
* неверные объяснения процессов. "обязательно убедитесь, что подключаете диод правильно <...> иначе толку от него будет ноль, и может сделать даже хуже." Человек вообще не понимает о чём пишет.

Найдите нормальный обучающий сайт, а лучше выбросьте эту ардуину в мусорное ведро.

Когда-то давно я занимался решением подобной задачи: хотел построить график зависимости силы притяжения от расстояния между двумя магнитами, магниты смоделировал просто как кольца с током.
1. Для вычисления силы использовал известную формулу, сила равна градиенту потенциальной энергии. Тут я сделал небольшое упрощение, рассматривал только силу, приложенную к центру магнитов - и не рассматривал то, что магниты могут пытаться повернуть друг друга. Но это тоже легко вычисляется.
2. Энергия магнитного поля, в свою очередь, это пространственный интеграл произведения B*H = mu*B*B по dV. Пространство ограничиваем какими-то разумными пределами, я использовал ~200-1000 диаметров кольца.
3. Магнитная индукция аддитивна, то есть B в "точке" dV складывается из B1 и B2 - полей одного и другого кольца. Кольца одинаковые, поэтому достаточно предварительно вычислить поле вокруг одного кольца, а потом для взятия интеграла складывать уже заранее вычисленные поля.
4. Магнитное поле кольца осесимметрично, поэтому достаточно вычислить поле всего лишь на одной четверти плоскости, построенной на радиусе кольца и его оси.
5. Для вычисления этого поля - формула Био-Савара-Лапласа, и циркулярный интеграл по контуру кольца. Тут уже наверное никак особо не соптимизируешь, разве что считать только половину кольца.

Таким образом, по сути задача сводится к: циркулярный интеграл Био-Савара-Лапласа по половине кольца, посчитанный в каждой точке четвертьплоскости, и интегрирование этого поля по объёму с нужными смещениями и поворотами.

В трёх местах есть проблемы:
1. Поле вблизи кольца. Есть хороший метод в учебниках физики (к сожалению, не помню где) - рассматривать не бесконечно тонкое кольцо, тогда эта проблема легко решается.
2. Грануляция сетки для более быстрого и точного расчёта. Стандартный способ - около кольца бьём плоскость на много сегментов, а дальше от кольца - сегменты всё больше и больше.
3. Несовпадение сеток двух колец, если кольца повёрнуты друг относительно друга - я бы попробовал алгоритм Брезенхема для поиска ближайшего сегмента или даже интерполяции. Повторюсь, повороты я не реализовывал.

Это было сделано лет пять назад на ноутбуке с Core2Duo 8500, на Delphi с применением SSE2, в 5 потоков. В зависимости от грануляции расчёт занимал до суток. В итоге была получена обратно-кубическая зависимость, как и ожидалось.

Если требуется определить поведение системы - устремляете систему к минимуму энергии вариацией положений/поворотов магнитов, я бы действовал методом Ньютона. Благо, пересчёт системы с новыми положениями не требует перерасчёта магнитного поля кольца, а лишь требует расчёта пространственного интеграла.

Что бы я сделал сейчас:
1. перенёс векторные вычисления (циркулярный интеграл) на CUDA, пространственный интеграл - наверное, тоже.
2. проверил, правомерна ли аппроксимация магнита кольцом с током - просто разместил бы много маленьких колец внутри объёма магнита (и конечно, в одном направлении, как домены), и сравнил поле такой имитации с полем одного кольца.
3. сделал полноценные графики градиентов энергии по всем степеням свободы. Думаю, нашлось бы много интересного.